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Abstract. In this contribution, we study the theoretical and numerical stabil¬ 
ity of a bidimensional relative velocity lattice Boltzmann scheme. These relative 
velocity schemes introduce a velocity field parameter called “relative velocity” func¬ 
tion of space and time. They generalize the d’Humieres multiple relaxation times 
scheme and the cascaded automaton. This contribution studies the stability of a 
four velocities scheme applied to a single linear advection equation according to 
the value of this relative velocity. We especially compare when it is equal to 0 
(multiple relaxation times scheme) or to the advection velocity (“cascaded like” 
scheme). The comparison is made in terms of L°° and L 1 stability. The L°° 
stability area is fully described in terms of relaxation parameters and advection 
velocity for the two choices of relative velocity. These results establish that no 
hierarchy of these two choices exists for the L°° notion. Instead, choosing the 
parameter equal to the advection velocity improves the numerical L 2 stability of 
the scheme. This choice cancels some dispersive terms and improve the numerical 
stability on a representative test case. We theoretically strengthen these results 
with a weighted L 2 notion of stability. 


Introduction 

Lattice Boltzmann schemes are applicable to many different fields such as hydrody¬ 
namics, acoustics, magnetohydrodynamics and multiphase fluids B EH E2 0 d, 
because the associated algorithm is simple, fast and flexible. This algorithm, associ¬ 
ated with a cartesian lattice and a finite set of velocities, consists in computing some 
particle distribution functions at discrete values of time: they are updated through 
two successive phases of transport (exact) and relaxation (also called collision). 

In the d’Humieres multiple relaxation times (MRT) framework [8], the relaxation is 
made thanks to a set of moments, linear combinations of the particle distributions. 
Each of these moments relaxes towards an equilibrium with a proper time scale. 
The MRT schemes have been widely studied in terms of consistency |10l RT| : 20 j but 
slightly in terms of stability [221 El- They can be particularly used to solve the 
Navier-Stokes equations [HI El SUSSES]. They however undergo some instabilities in 
the small viscosity limit [22]. 

A “cascaded automaton” [15], introduced in two and three dimensions for the Navier- 
Stokes equations, reduces these instabilities [16] . Some studies aim to understand 
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the algorithm of this automaton and its improvements in terms of stability. First, it 
has been written in the MRT framework with a generalized equilibrium depending on 
the non conserved quantities [Tj. It has then been included in a new class of relative 
velocity schemes mm that are the center of interest of the authors. Both MRT 
and cascaded are particular cases of the relative velocity schemes. 

These relative velocity schemes use a set of moments depending on a velocity held 
parameter called relative velocity for the relaxation. This parameter is a function of 
space and time. Their consistency has been studied for one and two conservation laws 
mm~ Their stability has also been investigated numerically in the case of the D 2 Q 9 
scheme for the Navier-Stokes equations m This study, using the L 2 von Neumann 
stability |221126] . shows the importance of the choice of the polynomials defining the 
moments, that justifies the stability gain obtained by the cascaded automaton. 

This contribution studies the theoretical and the numerical stability of a bidimen- 
sional scheme with four velocities for the simulation of a single linear advection 
equation. We expect to improve the stability for a relative velocity equal to the 
advection velocity (“cascaded like” scheme) instead of 0 (MRT scheme). This study 
is based on two stability notions. The first is a L°° notion that has been used for 
one dimensional and vectorial schemes [ 6 j 118] . The second is a weighted L 2 notion 
introduced in |3j. It has been applied to MRT schemes linearized around the zero 
velocity or in the BGK (Bhatnagar, Gross, Krook) case )4] for advection equations 
|3J EU [25]. 

In the first part, we briefly recall the framework of the relative velocity schemes. Then 
we present the four velocities scheme of interest. In the second part, the effect of the 
relative velocity on the stability is illustrated thanks to the method of the equivalent 
equations EH HO]. We particularly focus on the structure of the dispersion terms for 
two choices of relative velocity (0 or the advection velocity). These results are then 
illustrated by a numerical L 2 von Neumann stability study. The two last parts give 
theoretical stability results. In the third part, the L°° stability area is fully described 
in terms of relaxation parameters and advection velocity for the two choices of relative 
velocity. In the fourth one, a weighted L 2 notion (21j is used to validate the numerical 
results obtained in the second part according to the relative velocity. 

1. Framework 

This section presents first the relative velocity schemes in a general framework. It is 
then particularized to a bidimensional scheme with four velocities for the simulation 
of a linear advection equation. 

1.1. Presentation of the relative velocity schemes 

For d G N*, we consider a cartesian lattice £ of W d associated with a space step Ax. 
The space time is given by the acoustic scaling At = Ax/X for A £ M a velocity 
scale. A set of q € N* velocities denoted by v = {u 0 ,... depending on the 
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velocity scale is chosen such that for each node x E £, x + t>-A t is still a node of 
C. The lattice Boltzmann method computes the discrete values of several particle 
distributions denoted by Jf, 0 ^ j ^ q — 1 , the distribution f- being associated 
with the velocity v q . We denote by f = (f 0 ,..., f q _ 1 ) the vector of the particle 
distributions. An iteration of the algorithm consists in the succession of a phase of 
relaxation, nonlinear and local in space, and a phase of transport distributing the 
particles on the neighbouring nodes. 

The relaxation operator of the relative velocity schemes originates from the frame¬ 
work of the MRT schemes [ 8 j. Some moments, linear combinations of the particle 
distributions, relax with a proper relaxation rate. We define a matrix of moments 
depending on a velocity field u function of space and time. This matrix, supposed 
to be invertible, is characterized by a set of polynomials P k , 0 ^ k ^ q— 1, 

M{u) kj = P k {v ] - u), 0 ^ k,j ^ q- 1. 

Let us note that the MRT scheme corresponds to u = 0. The vector of moments is 
then obtained from the particle distributions through a linear transformation 

m(u) = M(u) f. (1) 

The relaxation then carries on the components m k (u), 0 ^ k ^ q — 1, of m(u ) = 
(m 0 (S),...,m q _ 1 (5)) 

ml(u) = m k (u)+s k (m e k q (u)-m k (u)), 0 < fc < q-1, (2) 

where m eq (u) are the moments at equilibrium and s k the relaxation parameters. Some 
relaxation parameters are chosen null to enforce some conservation laws. The equi¬ 
librium is supposed to derive from an equilibrium distribution / eq E R 9 independent 
of u 

m eq (u) = M(u)f eq . (3) 

The equilibrium f eq is an a priori non linear function of these conserved moments. 
The postcollision distributions are obtained by the linear transformation 

f* = M~ l (u)m*{u). (4) 

Then, the particles are advected from node to node 

fj{x, t + At) = f*(x - VjAt,t), 0 ^ j ^ q-1. 

In the following, the relative velocity u is a real constant. When u is specified, the 
relative velocity scheme is called scheme relative to u. 

1.2. The twisted relative velocity D 2 Q 1 scheme for a linear advection equa¬ 
tion 

In this section, we present a relative velocity D 2 Q 4 scheme for a linear advection 
equation. We expect the scheme relative to the advection velocity ( Tl = V), also 
called “cascaded like” scheme in the following, to be more stable than the MRT 
scheme u = 0 when this velocity V increases. 
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We want to approximate the following bidimensional advection equation 

d t p(x,t) + V ■ Xp(x,t) = 0, sGR 2 ,tGl + , (5) 

where V = ( V x , V y ) 6 M 2 is the advection velocity. With this aim, we consider 
the four velocity scheme called twisted D 2 Q 4 defined by the set of velocities v = 
{(A, A), (—A, A), (—A, —A), (A, —A)}, for A = Ax/ At the velocity scale (hgure[l]). The 
moments are associated with the polynomials 1, X, Y, XY. 


A 

'-S 



Figure 1 . Velocities of the twisted D2Q4 scheme. 


Another choice with four velocities could be the D2Q4 defined by the set v = 
{(A, 0), (0, A), (—A, 0), (0, —A)}, and by the moments associated with 1, X, Y, X 2 — Y 2 . 
However, we have chosen to work with the twisted scheme because this scheme is ex¬ 
actly a D 1 Q 2 2 in terms of velocities and moments. This feature simplifies the proofs 
of consistency and stability. Moreover, the stability results of the D2Q4 scheme can 
be deduced of the results for the twisted D2Q4 scheme: the stability areas correspond 
thanks to the composition of a rotation and a homothety (Appendix |A|) . 

The framework containing only one equation, we choose one conservation law on the 
density 

r = tfr 

3=0 

The other moments are relaxed with a relaxation parameter s ? gR for the first order 
moments and s xy € M for the second order one. To approach the equation (|5j), the 
equilibrium must read 

m eq (u) = ( p , (V x - u x )p, iyy - ^)p,m e 3 q (S)), ( 6 ) 

where m ^(u) is the second order moment equilibrium determining the diffusion 
terms. Two different equilibria are chosen: the non intrinsic equilibrium 

mf{u) = p(V x -u x )(V y -u y ), (7) 

and the intrinsic equilibrium 

m eq {u) = p(u x u y - u x V y - u y V x ). 


( 8 ) 
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The first equilibrium results from the will to see the twisted D2Q4 as the product 
of the D 1 Q 2 by himself: it could be seen as the product of two one dimensional 
equilibria. The second equilibrium is called intrinsic because it leads to a diffusion 
term that is independent of the basis of writing (proposition 2.2) but is not isotropic. 
It is important to note that these equilibria are chosen such that the relation (|3| is 
satisfied. Then the second order asymptotics do not depend on the relative velocity 

spang. 


2. Exploitation of the equivalent equations for the stability 

The purpose of this section is to predict the stability behaviour of the twisted scheme 
defined in the section [l] according to the choice of u. To do so, we present the third 
order asymptotics of the scheme: we discuss on the definite positivity of the diffusion 
tensor and on the dispersive third order terms. The discussion is then illustrated on 
a numerical study of L 2 stability. We show that the cascaded like scheme contrary 
to the MRT scheme cancels some dispersive terms and stabilizes the scheme in the 
same time. 


2.1. Discussion for the scheme with a non intrinsic diffusion 


We study here the third order equivalent equations of the twisted D2Q4 with a non 
intrinsic diffusion. The following proposition results from the particularization of a 
general derivation of the equivalent equations for a scheme with one conser¬ 

vation law, d,q € N* HU We are interested in the structure of the upper orders 
according to the velocity held u and their influence on the stability. 


Proposition 2.1 (Diffusion and dispersion operators). Given V e M 2 , ( s q , s xy ) e M 2 , 
the twisted D2Q4 scheme relative to a constant velocity u G M 2 associated with the 
equilibrium (’6|?) and with the relaxation parameters (0 ,s q ,s q ,s xy ), approaches the 
third order equation 


d t p{x. t) + V ■ Vp(x, t ) — A t V 2 : V 2 p(x, t ) + At 2 V 3 : V 3 p(x, t ) = (D(At 3 ), 

x<ER 2 ,teR + , (9) 

where the diffusion matrix T> 2 is given by 


V 2 = a q 


(A 2 - {V x ) 2 ) 0 

0 (A 2 -^) 2 ^’ 


( 10 ) 


and the dispersion matrix V> 3 by 

A 2 - m 2 )(l - 12a 2 ) 2a q (X 2 —(V y ) 2 )(u x —V x )(a q —c 
2a q {\ 2 -{V x ) 2 )(uy-V y ){a q -a xy ) ^(A 2 - (P^) 2 )(l - 12a 2 ) 


V 3 = 


• ( 11 ) 
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The Henon’s parameters reas a q = l/s q — 1/2 and a xy = 1 /s xy — 1/2, the second and 
third order derivative operators are given by 

V 2 p= ($p v 3 P =(f p ^ y 3 yp ), 

\d£ y p cfip) \8§xyP 9yP J 

and : is the scalar product for the matrices viewed as vectors o/M 4 . 



These equations are obtained thanks to a formal calculus software, that implements 
an algorithm to obtain the equivalent equations in the linear case |2]. The following 
lemma exhibits the area of well-posedness of the second order truncation of ©• 


Lemma 2.1 (Well-posedness of the second order equation). The diffusion matrix V 


defined by (10) belongs to the space S. 
if and only if IV^ < A. 


++ 


of symmetric positive definite matrices 


Proof. The spectrum of the matrix V> 2 is {^ 2 — (V X ) 2 ,X 2 — ( V y ) 2 } that closes the 
proof. □ 


Because of the relation (J3|, the velocity parameter u appears only at the third order 
of the equations. Let us discuss of the structure of these third order dispersive terms 
according to the choice of u. 

First we choose the velocity field u equal to the advection velocity V (“cascaded 
like” scheme). The dispersion terms are then small compared to the diffusion terms. 
Indeed, this choice cancels the two off-diagonal components of T> :>) depending on the 
parameter a xy . Thus resting dispersion terms are W(A 2 — (V x ) 2 )(l — 12cr 2 )/6 and 
its symmetric in y. When a q tends to 0, those terms are equivalent to V (A 2 — 
{V x ) 2 ) multiplied by a constant. Thus when the diffusion (A 2 — (V x ) 2 ) decreases, the 
dispersion decreases at least with the same speed. 

On the contrary, the dispersion terms of the MRT scheme (u = 0) create instabilities 
when the diffusion (represented by cr q ) is weak. Indeed, the off-diagonal terms are 
conserved for this choice of velocity parameter: they are given by —2cr q V y (X 2 — 
( V x ) 2 )((J q — <J xy ) and its symmetric in y. This term is equivalent to 2a xy a q V y (X 2 — 
( V x ) 2 ), when a q gets close to 0. The dispersion terms are thus depending on the size 
of the parameter a . If a xy is important and o q tends to 0, the numerical stability 
should be deteriorated by dispersion phenomena for V x or V y close to A: indeed, 
the dispersion, behaving as 2a xy (j q X(X 2 — (L x ) 2 ), becomes greater than the diffusion 
given by a ( A 2 — ( V x ) 2 ). Instead, taking a xy close to 0 should limit the dispersion 
effects. 


Remark 2.1 (Exact third order scheme). We note that the choices u = V and 
(T q = 1/Vl2 cancel the dispersion matrix T> 3 given by (11). Thus the scheme is 
consistent with the second order advection diffusion truncation of 0 at the third 
order. Moreover, there is still a degree of freedom with the parameter a xy . This is 
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impossible for the MRT scheme (u = 0) excepted when there is only one relaxation 
parameter (BGK scheme: a q = o xy ). 

2.2. Discussion for the scheme with an intrinsic diffusion 

We discuss on the structure of the third order equivalent equations for the twisted 
D 2 Q 4 intrinsic scheme obtained thanks to the linear algorithm derived in |2j. The 
study is completely analogous to the one of the section [2~Tj 


Proposition 2.2 (Diffusion and dispersion operators). Given V £ M , (s q , s xy ) £ 
M 2 , the twisted D 2 Q 4 scheme relative to a constant velocity u associated with the 
equilibrium (7m8) and with the relaxation parameters ( 0 , s q , s q , s xy ), approaches the 
ra 


following third order equation 

d t p(x, t) + V ■ S7p(x, t ) — A t V 2 : X/ 2 p(x, t ) + At 2 V 3 : X7 3 p(x , t ) = 0(At 3 ), 


X £ 


1 1 £ 


where the diffusion matrix T> 2 is given by 

= ( <Jq (\ 2 — (V x ) 2 ) —G q V X V V 

2 V -o q v x vy a q {\ 2 - (vy) 2 ) J ’ 

and the dispersion matrix V 3 by 

(fAv x ) cj) 2 (u x ,u« ,v x ,vy) 

Mw,u x ,vy,v x ) Mv y ) 




(13) 


V 3 = 


with 


yx 


4 > i(.V x ) = —(A — (y x Y)(l — 12a q ), 


<j) 2 (u x ,u y , V y ) = - ( - V x (V y ) 2 + 4u 2 (A 2 {u x - V x ) + 2>V x (y y ) 2 - u y V x V y 

-u x {V y ) 2 ) - 4a q a xy (X 2 (u x - V x ) - (V x u y V y + d x {V y ) 2 ). 

The Henon’s parameters a q and a ar e given by 1 /s q — 1/2 and 1 /s xy — 1/2, the 
operators V 2 and V 3 are given by (12) and : is the scalar product for the matrices 
viewed as vectors 0 /M 4 . 

Remark 2.2. As expected, the diffusion is intrinsic. It is independent of the basis 
of writing. Indeed it is equal to 

^2 : V 2 = X 2 (d x + d 2 ) — (V ■ V) 2 . 

Lemma 2.2 (Well-posedness of the second order equation). The diffusion matrix T> 2 
defined by (13) belongs to the space 5^" + (M) of symmetric positive definite matrices 
if and only if \V I2 < 


Proof. The eigenvalues of the matrix T> 2 are ^ an( i — |V| 2 , that closes the proof. 

□ 
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Figure 2. Well posedness areas in V for the twisted D 2 Q 4 schemes: 
with intrinsic diffusion (circle of radius A), with a non intrinsic diffu¬ 
sion (square [—A, A] 2 ). 

Let us compare the well-posedness areas of the schemes with a non intrinsic diffusion 
(section |2.1| ) and an intrinsic diffusion. These areas are represented on the figure [ 2 ] 
We note that the scheme with a non intrinsic diffusion has a greater area in V. 
Consequently, it should authorize larger stable velocities than the scheme with an 
intrinsic diffusion. 

We now get back to the intrinsic case and discuss on the dispersion terms. As in the 
non intrinsic case, the scheme should be stable for larger velocities with the choice 
u = V than with u = 0 . Indeed, when a q tends to 0, the term </> 2 (u x , u y , V x , V y ) 
is of third order in V if u = V. Instead, if u = 0 , it is of first order in V. More 
dispersion is created for u = 0 when V x or V y are close to A and a xy is big. 

2.3. Illustration with numerical stability experiments 

In this section, we study the numerical L 2 stability for the twisted D2Q4 scheme 
according to the choice of the velocity parameter u. The link with the previous 
consistency study is made: the third order dispersion terms have an influence on the 
numerical stability. The choice u = V provides larger stability areas in V. 

We present a first numerical experiment with the advection of a circular spot: the 
initial conditions are given by 

p(x, t) = 1 + l c (®, t), x e [0, l] 2 , t 6 M, 

where C is the disc centered in (1/2,1/2) of radius 0.1 and l c is the function equal 
to 1 on this disc, null elsewhere. This spot moves with an advection velocity kGl 2 
in the domain [0, l] 2 constituted of 128 2 points with periodic boundary conditions. 
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a 1 

1/10 

1/20 

1/50 

1/100 

1/200 

D 2 Q 4 non intrinsic u = 0 

1.00 

0.80 

0.49 

0.34 

0.23 

D 2 Q 4 non intrinsic u = V 

1.00 

1.00 

1.00 

1.00 

1.00 

D 2 Q 4 intrinsic u = 0 

1.00 

0.79 

0.48 

0.33 

0.23 

D 2 Q 4 intrinsic u = V 

1.00 

1.00 

1.00 

1.00 

1.00 


Table 1. Maximal |V| stable in A scale for 0 = 0, a xy = l/\/3 and 
128 2 points. 



1/10 

1/20 

1/50 

1/100 

1/200 

D 2 Q 4 non intrinsic u = 0 

1.41 

0.80 

0.42 

0.28 

0.20 

D 2 Q 4 non intrinsic u = V 

1.41 

1.41 

1.41 

1.41 

1.41 

D 2 Q 4 intrinsic u = 0 

0.75 

0.56 

0.36 

0.26 

0.18 

D 2 Q 4 intrinsic u = V 

0.86 

0.76 

0.65 

0.59 

0.53 

Table 2. Maximal 

stable in 

A scale for 0 

II 

_£>• 

= 1/\/3 and 


128 2 points. 


The scheme is considered as stable if it has not broken after 2000 iterations. We 
present the biggest advection velocities keeping the scheme stable for different choices 
of relaxation parameters. Two different directions of advection velocities are consid¬ 
ered: 0 = 0 and 0 = 7 t/ 4 . We choose different values for the parameter a q . the 

parameter a xy being set to 1/ V3. The results are presented in the tables [l] and [ij 

The scheme relative to the advection velocity u = V authorizes larger stable ve¬ 
locities than the MRT scheme (u = 0 ). Most of the time, its stability areas do not 
depend on the relaxation parameters while they are bounded by 0 and 2. Instead, the 
stability areas of the MRT scheme decreases when cr q tends to 0. This phenomenon 
occurs for the scheme relative to u = V with an intrinsic diffusion in the direction 
0 = vr/4: however the instabilities appear for larger V than the MRT scheme. 

The dispersion terms exhibited in the equivalent equations (section [2j originate these 
phenomena. Their presence causes instabilities when a q is too small compared to a xy 
and | V|. We have seen that these terms cancel for u = V but are present when u = 0 . 
This explains the Constance of the stability area for u = V and its deterioration for 
u = 0 when a q decreases. 

These dispersive phenomena are illustrated on the figures [3] and [4] for the scheme 
with an intrinsic diffusion. We choose a xy = l/\/3, cr q = 1/20 and V = (0.9A,0) for 
these draws. According to the data of the table [lj this configuration is a stable case 
for the scheme relative to u = V and an unstable case for the MRT scheme. The 
spot is represented at the times t = 0 and t = 0.4 for 256 2 points. 
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Figure 3. Advected spot of velocity V = (0.9A,0) for the twisted 
D 2 Q 4 scheme relative to u = 0 (MRT) with an intrinsic diffusion at 
t = 0 (left) and t = 0.4 (right) for a xy = l/\/3, cr q = 1/20. 



Figure 4. Advected spot of velocity V = (0.9A, 0) for the twisted 
D 2 Q 4 scheme relative to u = V with an intrinsic diffusion at t = 0 
(left) and t = 0.4 (right) for a xy = 1/Vs, a q = 1/20. 


The dispersion for the MRT scheme is clearly visible on the figure [3] some oscillations 
appear behind the advected spot. The scheme is going to break since the density p is 
increasing. These phenomena are absent for the scheme relative to u = V (figure [4]) 
that remains numerically stable. 


Our second numerical experiment uses the L 2 linear stability of von Neumann. We 
discuss on the spectrum of the amplification matrix of the scheme. This matrix, 
characterizing an iteration of the scheme, has to be first determined. The equilibrium 

and 


being linear (6|7|8 1, there exists a matrix E = E(V,s ) for s = (0, s c 


’qi Sqi &xy) 


V £ 


so that 


f eq = Ef. 

The relaxation phase of the relative velocity D 2 Q 4 scheme reads 
/* = (/ + M(u)~ 1 DM(u)(E - I))f , 


where D = diag(s) is the diagonal matrix of the relaxation parameters. This expres¬ 
sion holds for each node x of the lattice, the relaxation being local in space. The 
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Figure 5. Velocities V stable in A scale for the scheme relative to 
u = 0 (MRT on the left), u = V (on the right), for a non intrinsic 
diffusion with s n = 1 and = 1. 



expression of the distribution after the transport phase is given by 
f j (x,t + At) = [{I + M(u)~ 1 DM(u)(E-I))f} j (x-v j At,t), xe£,teR. (14) 

Taking the Fourier transform of ( |14| ), the transport operator becomes local in space 
and is represented by the diagonal matrix A = A(k) for k 6 M 2 whose diagonal 
components are given by e l ^ tk - v j ; 0 ^ j ^ 3. The amplification matrix then reads 
L(u) = L(u,V,k,s ) = A(I + M(u)~ 1 DM(u)(E — /)). It characterizes a time 
iteration of the scheme in the Fourier space 

f(k,t + At) = L(u)f(k,t), teR, 

where f is the Fourier transform of f. We want to exhibit the advection velocities 
V for which the scheme verifies the necessary condition of L 2 stability 

max r(L(u )) ^ 1, (15) 

fceM 2 

where r is the spectral radius. We are thus interested in the set 

{V = (V x , V y ), max r(L(u)) < 1}. 
fees 2 

The figures [5] to [lO] present this set in the plan (V x , V y ) for the twisted D 2 Q 4 scheme 
with a non intrinsic diffusion and some relaxation parameters s q and s xy . The left 
draw is about the MRT scheme (u = 0 ) and the right one illustrates the case u = V. 
We expect the right areas to be larger than the left ones. 

Firstly, the draws corresponding to a BGK scheme (s q = s xy = 1) are independent of 
u (figure [5j). This result was expected because the velocity held u does not appear 
in the scheme for one single relaxation parameter since ([3]) is verified. Secondly, the 
scheme relative to u = V verifies the necessary condition of stability on the biggest 
area \V 1^ ^ A (CFL) for all the s chosen (figure [ 5 ] to [To] ). The Constance and 
the optimality of these areas confirm the phenomena observed for the advected spot 
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Figure 6 . Velocities V stable in 
u = 0 (MRT on the left), u = V 
diffusion with s n = 1 and = 1.5. 
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Figure 7. Velocities V stable in A scale for the scheme relative to 
u = 0 (MRT on the left), u = V (on the right), for a non intrinsic 
diffusion with s„ = 1 and s„„ = 1.9. 

h 


(tables [I] and [2]). Finally, the stability areas of the MRT scheme decrease as the two 
relaxation parameters are moving away from each other. The stability area is reduced 
to V = 0 for s n = 2 or = 2. Whatever the choice of s, the areas of the MRT 
scheme are included in those of the scheme relative to u = V. These results are also 
consistent with the test case of the advected spot: the deterioration of the areas are 
due to the third order dispersive terms of the equivalent equations (proposition | 2 . 1 [). 


For the D 2 Q 4 scheme with an intrinsic diffusion (figures 11 to [16| ), the results are 
different but the general trend is the same. The scheme relative to u = V verifies 
(15) for larger sets of V than the MRT scheme, excepted when s g = 1, s xy = 1.5 


(figure 12). However, the areas associated to u = V are not optimal any more: they 
decrease when s q and s xy move away from each other. Both velocity fields (u = 0 and 
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Figure 8 . Velocities V stable in A scale for the scheme relative to 
5 = 0 (MRT on the left), 5 = V (on the right), for a non intrinsic 
diffusion with s„ = 1 and = 2. 


0.8 

0.6 

0.4 

0.2 

0 


-0.2 


-0.4 


-0.6 


-0.8 


—1*— 
-1 


T 

. ... 0.8'- 

. ... 0.6'- 

. ...... 0.4 - 

.. 0.2 

.;;;;;;;;;;;. > °i .............. . 

.. -0.2- 

. . -0.4- • ........................................................................ 

. ... -0.6- 

- 0 . 8 - "»- - * " — 

-i-i-1-i -li—.—.—.—.— i —.—.—.—.—A—.—.—.—.—i—.—.—.—.—A 

-0.5 0 0.5 1 -1 -0.5 0 0.5 1 

VX vx 

Figure 9. Velocities V stable in A scale for the scheme relative to 
5 = 0 (MRT on the left), 5 = V (on the right), for a non intrinsic 
diffusion with s n = 1.9 and = 1. 

h 



u = V ) undergo this phenomenon but the effect of the dispersion is bigger for the 
MRT scheme. These results confirm the advected spot ones: the same deterioration 
appears when a q tends to 0 (table |TJ and |2j) to a lesser extent for u = V than for 

u = 0 . 

The main conclusions of this section are the following. Choosing u = V cancels some 
dispersive terms becoming important as V grows. This cancellation does not occur 
for 5 = 0 . We must view these results in parallel with the fact that the scheme is 
more stable (numerical L 2 notion) for u = V than for 5 = 0 . 
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Figure 10. Velocities V stable in A scale for the scheme relative to 

u = 0 (MRT on the left), u = V (on the right), for a non intrinsic 

diffusion with s n = 2 and = 1 . 
h 


vx vx 

Figure 11. Velocities V stable in A scale for the scheme relative to 
u = 0 (MRT on the left), u = V (on the right), for an intrinsic 
diffusion with s n = 1 and s .= 1 . 

h 


3. Theoretical L°° stability 

We study in this section the stability of the relative velocity D 2 Q 4 scheme with 
respect to a notion of L°° stability presented in [ 6 , 18]. We use it to demonstrate 
some maximum principles for u equal to 0 (MRT scheme) and V, the advection 
velocity (“cascaded like” scheme). The L°° stability area is fully described in terms of 
relaxation parameters and advection velocity for these two choices of relative velocity. 
We show that this L°° notion does not differentiate u = 0 from u = V in terms of 
stable behaviour. In each case, the parameters (s , s xy ) corresponding to a non empty 
area of L°° stability in V are included in the square [0, 2] 2 . 
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Figure 12. Velocities V stable in A scale for the scheme relative to 
u = 0 (MRT on the left), u = V (on the right), for an intrinsic 
diffusion with s n = 1 and = 1.5. 
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Figure 13. Velocities V stable in A scale for the scheme relative to 
u = 0 (MRT on the left), u = V (on the right), for an intrinsic 
diffusion with s„ = 1 and = 1.9. 

h ^y 


Definition 3.1 (L°° stability). Let us consider a D c {Q q scheme on the cartesian 
lattice C, £ M d . The mass at time t is given by 

q-l 

p to \t) = ^2fj (*>*)> teR - 

x£C j =0 

Let C £ R* , suppose that the initial distributions of particles are nonnegative and 
that the initial mass is bounded by C: 

fj(x, 0)^0, 0 ^ j ^ q - 1, x £ £, and p tot (0) ^ C, 

the scheme is said to be L°° stable if for all time t n £ R, n £ N, 

f-(x,t n ) ^0, 0 ^ j ^ q — 1, x £ £, and p tot (t n ) ^ C. 
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Figure 14. Velocities V stable in A scale for the scheme relative to 

u = 0 (MRT on the left), u = V (on the right), for an intrinsic 

diffusion with s n = 1 and = 2. 
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Figure 15. Velocities V stable in A scale for the scheme relative to 
u = 0 (MRT on the left), u = V (on the right), for an intrinsic 
diffusion with s„ = 1.9 and = 1. 

H a '£/ 


The proofs of this a priori non linear L°° stability notion are only depending on the 
relaxation phase. Indeed, the transport exchanges the distributions between each 
other: it does not act on their positivity and on their mass. Thus it is sufficient to 
show that the postcollision distributions remain nonnegative and the mass bounded. 
Note also that if p tot (0) is bounded, p tot (t n ) remains bounded for all time t n because 
this mass is conserved by the scheme. So there is only to study the positivity of the 
postcollision distributions functions. 

To do so, we express the postcollision distributions of particles /* as a function of 
the precollision one f for a general D c iQ q scheme with one conservation law. These 
postcollision distributions are first expressed as functions of the postcollision moments 
thanks to Q. After applying the relaxation identity Q, the precollision moments 
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Figure 16. Velocities V stable in A scale for the scheme relative to 
u = 0 (MRT on the left), u = V (on the right), for an intrinsic 
diffusion with s n = 2 and = 1 . 


are written in the frame of the particle distributions with ([Tj) . These calculations 
lead to the following expression of the relaxation. 

9—1 9-1 9-1 

fi* = (i - 2 &%**&)* M &)k 3 )fj 

k =1 j =0 k =1 

9-1 

+ s k M (uXk m T&), 0 < l < q - 1 . (16) 

k =l 


The framework of study is linear because we are approaching a linear advection 
equation. As a consequence, the equilibrium term of (16) is a linear function of 
the conserved variable p and thus of the distributions f-. Consequently, there is a 
matricial relation between /* and f : sufficient conditions of L°° stability are obtained 
checking when this matrix is nonnegative. 


Definition 3.2. For q 6 N*, a square matrix of is said to be nonnegative if 

all its coefficients are nonnegative. 


Note that this notion is different from a positive matrix that have all its eigenvalues 
nonnegative [23]. In the following, we prove the results in the twisted case: the 
analogous ones for the D 2 Q 4 scheme are available in the appendix [B] and use the 
proposition A.l of the appendix |A| We first consider the relative velocity D2Q4 scheme 
with a non intrinsic diffusion (equilibrium ([7])). The proofs are quite technical and 
involves geometric inequalities on V and s. 
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3.1. The non intrinsic diffusion case 


Proposition 3.1 (L°° stability areas for the MRT scheme). Let V S M , (s q , s xy ) E 
10)2 consider the twisted D 2 Q 4 MRT scheme (u = 0) with the relaxation parameters 

1 ? 


( 0 , s q , s Q , s xv ), associated with the equilibrium 


m eq (0) = p(l,V x ,V y ,V x V y ). 


Note V = (V' T , V y ) where V x = V x + V y , V y = V x - V y and 7 = s q /s xy : 

• if 0 < s xy ^ min(s q ,2 — s q ) (area BCD on the figure [7?j) ; the scheme is L°° 


stable for all V so that 

(V' T ± 2A 7 ) 2 - (V y f ^ 4A 2 ( 7 2 - 1), (17a) 

(V y ± 2 A 7) 2 - (V *) 2 ^ 4 A 2 ( 7 2 - 1). (17b) 

if s q ^ s xy ^ 2 s q and s q ^ 1 (area ABC on the figure |7t|), the scheme is L°° 
stable for all V so that 

{V x ± 2A 7 ) 2 — (V y ) 2 ^ 4A 2 ( 7 — l) 2 , (18a) 

(V y ± 2 A 7) 2 — (V x ) 2 ^ 4A 2 ( 7 — l) 2 . (18b) 

if2 — s q ^ s xy ^ 2(2 — s q ) and s q ^ 1 (area ACD on the figure [7t[) ; the scheme 
is L°° stable for all V so that 

( V x ± 2A 7 ) 2 - (V ^) 2 ^ 4A 2 ((7 + l ) 2 - —) , (19a) 

v s xy J 

{V y ± 2A 7 ) 2 - (V x ) 2 ^ 4A 2 ((7 + l ) 2 - —). (19b) 

xy 


• if s xy = 0 and 0 < s q ^ 2 (ray ]BD] on the figure |l?l), the scheme is L°° stable 
for V = 0 . 

• */ s xy = s q = 0 (point B on the figure |~/7p, the scheme is unconditionally L°° 
stable. 

• For all other s, there is no V corresponding to a L°° stable scheme. 


In particular, the parameters (s q , s xy ) corresponding to a non empty area of L°° sta¬ 
bility in V are inclu,ded in the square [0, 2] 2 . 


conditions. 


17 


The limit 


The stability areas in the plan (s , s xy ) are represented on the figure 
cases for the parameters s and s lead to consistent inequalities for' the velocity 


Proof. The first step consists to express /* as a function of f thanks to the relation 
(16). If the matrix sending f on f* is nonnegative, the scheme is L°° stable. A formal 
calculus software guarantees that it corresponds to ensure the following relations: 

/ S xv 2 s a — s rv 2 s + s„ 

min ( — —--^,1-? 

v 4 4 


'xy 


)±^{v x +(-iyv y )+(-iy^v x v y >o, * = o,i. 


4A 2 


4 


4 
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Figure 17. L°° stability areas in s of the twisted D 2 Q 4 MRT scheme 
with anistropic diffusion. BCD : (17a|17b), ABC : (18a|18b), ACD : 
( |19a|19bp . 


If s xy = 0, these relations imply V = 0 and 0 ^ s q ^ 2. Otherwise, we factorise by 
s xy /A A 2 and express these identities as functions of the coefficient 7 = s q /s xy . This 
leads to 

A 2 min (1,27-1,——27-I) ±X-f(V x +(-iyV y ) + {-iyV x V y ^ 0, i = 0,l, (20) 

s xy 

if s xy > 0 or 

A 2 max (1, 27 - 1, — - 27 - 1) ± X^(V X + (-1 yv y ) + (-1 fV x V y < 0, * = 0,1, 

s xy 


if s xy < 0. These inequalities are expressed thanks to quadratic forms in V. To 


characterize the associated geometry, we write them in an adapted basis. To do so, 
we use the variable change 

( V 2 ') 2 - ( V ^) 2 

Y x V y = - 1 —— - v ’ for V x = V X + V y , V y = V X -V y . 


The inequalities ( |20[ ) then read 

4A 2 min (l, 2 7 - 1, — - 2 7 - l) ± 4A7V 1 + (V x ) 2 - (y y f ^ 0, 


xy 


completed by the ones obtained when the roles of V x and V v are exchanged. We then 
center them in the adapted frame to obtain if s xy > 0 
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FIGURE 18. Value of the maximum of (21). Area 1 
(7 + l ) 2 - 4/s area 3 : (7 ~ l) 2 - 


1 , area 2 : 


(V x ± 2 A 7) 2 - (y y ) 2 ^ 4A 2 max (7 s - 1 , (7 - l) 2 , (7 + l ) 2 - —), ( 21 ) 

s xy 

plus the inequalities obtained exchanging the roles of V x and V y . The case s xy < 0 
yields to 

(V* ± 2A 7 ) 2 - {V y f < 4A 2 min ( 7 2 - 1 , (7 - l) 2 , (7 + l ) 2 - —). ( 22 ) 

S'r'u 


We begin by showing that there is no velocity stable in the case s xy < 0. The 
minimum of ( 22 ) is equal to (7 — l ) 2 (nonnegative) if s q ^ s xy < 0 , to 7 2 — 1 (negative) 
if s xy ^ min(s , 2 — s q ) and (7 + l ) 2 — 4/s (nonnegative) if 2 — s ^ s xy < 0. In 


the first and third cases, (| 22 j) would have a solution if and only if the left pole of the 
hyperbole centered in 2 A| 7 | is negative. We then should have 

2A|y| — 2A|7 — 1| ^ 0 , 


in the first case and 

2A| 7 | — 2 A (( 7 + l ) 2 — 4/ s xy )^ ^ 0 , 

in the third one. Considering that s xy < 0, the first case is equivalent to s xy ^ 2 s q 
that has no intersection with s q ^ s xy < 0 where (7 — l ) 2 is the minimum. The third 
one is equivalent to s xy ^ 2(2 — s q ) whose intersection with 2 — s q ^ s xy < 0 is empty. 
It remains the case 

(yyf _ (y* ± 2 A 7) 2 ^ 4A 2 (1 - 7 2 ), 

when s xy ^ min(s (J ,2 — s ). These inequations combined with the ones obtained 
exchanging the roles of V x and V y have no common intersection. This closes the case 

®xy V 
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FIGURE 19. Layout of the bounds of the stability areas in V for the 
twisted D2Q4 MRT scheme. 


We exhibit now the conditions for which (21) has at least a solution V. The eventual 


stability area is located between the two hyperboles 

(V x ± 2Ay ) 2 - (V y ) 2 = 4A 2 max (y 2 — 1, (7 — l) 2 , (7 + l) 2 - — ), 

s xy 

plus the ones obtained exchanging the roles of V x and V y . Let us note that this maxi¬ 
mum is always nonnegative because (7— l) 2 is nonnegative. This area is delimited by 
the points A, B , C, D on the figure |l9} The inequalities (21) have a solution V if and 
only if the left pole P of the hyperbole centered in 2 A 7 has a nonnegative abscissa 
(figure 19). When the maximum is y 2 — 1, this abscissa is 2A(y — I 7 2 — 112 ) 5 that is 


nonnegative if s q ^ 0. When it is (7 — 1) , the abscissa 2 A (27 — 1) is nonnegative if 
and only if s xy SC 2 s q . Finally when the maximum is (7 + l ) 2 — 4/s,.,,, the abscissa is 


xy > 


equal to 2 A (7 — ((7 + 1)“ — 4/s xy ) 2 ). The positivity of this coefficient is equivalent 
to s xy < 2(2 - sl). 


It remains to determine the maximum of (21) in the case where there is at least one 


velocity V stable: three cases, represented on the figure 


18 


appear. 


If 5. 


'xy 




2 — s xy , the maximum is 7 — 1 , if a. 
s q < mm(s xy ,l), it is ( 7 - l) 2 . 


'xy 


^ 2 - s q and s q ^ 1 , it is (7 + 1 )“ - 4 /s xy , if 


□ 


From the proposition 3.1 we deduce the simpler case of the BGK scheme correspond¬ 


ing to one single relaxation parameter. 
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Corollary 3.1 (The BGK case). Let V £ M 2 , u £ M 2 , s £ M, consider the twisted 
D 2 Q 4 scheme relative to u, BGK of parameter s, associated with the equilibrium (6^1 

m e *{u) = p( 1, V x - u x , V y - u y , ( V x - u x )(V y - u y )). 

Note V = (V x , V y ) where V x = V x + V y , V y = V x - V y : 


• if 0 ^ s ^ 1, the scheme is L°° stable for IV^ ^ A. 

• if 1 ^ s ^ 4/3, the scheme is L°° stable for all V so that 

(V x ±2A) 2 -(V s 0 2 ^16A 2 (l- (23a) 

(V y ±2A) 2 -(V x ) 2 ^16A 2 (l-i). (23b) 

• If s > 4/3, no V corresponds to a L°° stable scheme. 


In particidar, the parameters s corresponding to a non empty area of L°° stability in 
V are included in [ 0 , 2 ]. 

Remark 3.1. The BGK scheme does not depend on the velocity field u since the 
relation 0 is verified. 


Proof. In this context, the parameter 7 is equal to 1. 
summarize in 


Thus the inequalities (21) 


{V x ± 2A ) 2 - ( V y ) 2 > max (0,16A 2 (l - J)), 

plus the inequality obtained exchanging the roles of V' r and V y . If s ^ 1, it is 
equivalent to 


(V x ± 2A ) 2 - {V y ) 2 ^ 0, 


and so to IVI^ ^ A, else we obtain (23a p3b| . 


The inequalities ( 23a|23b ) admit a non empty set of stable velocities V if the abscissa 
of the left pole of the hyperbole 


(V x - 2A ) 2 - ( V y ) 2 = 16A 2 (l - i) , 

is nonnegative. This abscissa, equal to 2A(1 — 2(1 — 1 /s) 2 ) for s ^ 1, is nonnegative 
if and only if s ^ 4/3. □ 


For s ^ 1, the BGK scheme is optimal in terms of L°° stability (CFL reached). When 
s q and s xy go far from each other (TRT scheme), the stability area decreases until it 
reaches V = 0 at the boundaries of the triangle represented on the figure 17 There 
is no stable velocities outside of this triangle. 


We now focus on the scheme relative to u = V with a non intrinsic diffusion and two 
relaxation parameters. We show that the optimal area of stability spreads on some 
TRT schemes contrary to the MRT scheme. 
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Proposition 3.2 (L°° stability areas for a relative velocity scheme). Let V E 


0, 


'qi ® xy) ^ 


consider the twisted D 2 Q 4 scheme relative to u = V associated with 


the relaxation parameters ( 0 , s ( 


'qi Sq> $xy) ’ 


and 


the equilibrium 03 


m e yV) = p(l, 0,0,0). 


Let us note V = (V x , V y ) where V x = V x +V y , V y = V x — V y and^y = s xy /(2s q — s xy ), 
then: 


if s q ^ s xy ^ min(l, 2 s q ) (area ABC on the figure^2(fy, the scheme is L°° stable 


on 


^ X. (24) 

if s xy < 2 s q , max(s 5 ,1) ^ s xy ^ 2(2 — s ) (area BCED on the figure 20), the 
is L°° stable on the intt 

(' V x ± 2A 7 ) 2 - (V y ) 2 > 


scheme is L°° stable on the intersection of (24) with 

16A 2 


(V y ± 2A 7 ) 2 - {V X Y > 


(2 S q S X y) 

16A 2 
(2s 0 - a 


(a*y s g (2 Sg)), 

2^xy ®q(2 Sq)). 


(25a) 

(25b) 


xy> 


if 0 ^ s xy ^ min(Sq, Sq (2 — s )) farea ACF on i/ie figure 20), the scheme is 


xy 

L°° stable on 


< Ay. 


(26) 


if s q (2 — s q ) ^ ^ min(sg, 2(2 — s )) (area CEF on the figure 20), the scheme 

is L°° stable on (25a\25b ). 

For all other s, no V corresponds to a L°° stable scheme. 

if s xy = 2s q and 0 < s xy ^ 2 (ray ]AD] on the figure 20), the scheme is L°° 

stable on the intersection of p 4 | ) and 

'2 - s„ 


\V\ 1 sf 2A 


xy 


xy 


(27) 


• if s xy = Sq = 0 (A on the figure \2(ty , the scheme is unconditionally L°° stable. 

In particular, the parameters (s q ,s xy ) corresponding to a non empty area of L°° sta¬ 
bility in V are included in the square [0, 2] 2 . 


Proof. We determine the sufficient conditions for the positivity of the matrix sending 
f on /*. A formal calculus software leads to the following inequalities 


(2Sq — s xy )(X ± V X )(X ± V y ) 0, (28) 

(A ± V x )(±(2s q - s xy )V y + Xs xy ) 0, (29) 

AX 2 -((2 s q + s x y )X 2 ± s X yX(V x +(-iyV y )+(-l)fi s xy-2sq)V x V y ) fiO, i = 0,1, (30) 


plus the inequality coming from ( |29| ) when the roles of V x and V y are exchanged. 
We begin by the study of (28). 
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Figure 
to u = 
BCED : 


20. L°° stability area in s of the twisted D 2 Q 4 scheme relative 
V with a non intrinsic diffusion. ACF : (26), ABC : (24), 
fl25a|25bp , CEF : fl25a) . 


Case 1 : 2 s q < s xy . These inequalities are equivalent to 

(. \±V x )(X±V y ) < 0. 

No V verifies this inequality: such a set of V would be the intersection of fourth 
areas that do not intersect. We must assume that s xy ^ 2 s q to eventually have a non 
empty stability area in V. 


Case 2 


Sxy 


If s xy = 0 , all the inequalities are obviously verified and the scheme is unconditionally 
stable in V. If s xy < 0, A ± V x ^ 0 that is impossible. If s xy > 0, (29) becomes 
(A ± C' T )A ^ 0, plus the analogous inequality in V y : this is equivalent to (24). The 
inequalities (30) read 


1- -^-{2X±(V X ±V y )) ^0, 
4A 


that is equivalent to (27). The final stability area is the intersection of (24) and (27). 
Three cases, represented on the figures [2lJ [22] and 23, are then possible. 


Let s xy ^ 1, then 


2A 


s xy ) 


> 2A, 


x y 
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Figure 21. L°° stability area for the scheme relative to u = V with 
a non intrinsic diffusion for s xy = 2 s q . The point P is of abscissa 
2A(2 — s xy )js xy ^ 2A. 



Figure 22. L°° stability area for the scheme relative to u = V with 
a non intrinsic diffusion for s xy = 2 s q . The point P is of abscissa 
A < 2A(2 - s xy )/s xy ^ 2A. 


and the stability area is given by (24) (figure 21), let 1 s xy ^ 4/3, then 


(2 s x 


A ^ 2A--^ ^ 2A, 

s xy 

















26 


FRANQOIS DUBOIS, TONY FEVRIER, AND BENJAMIN GRAILLE 



Figure 23. L°° stability area for the scheme relative to u = V with 
a non intrinsic diffusion for s xy = 2s . The point P is of abscissa 
2A(2 — s xy )/s xy A. 


22 


and the stability area is the intersection of (24) and (27) represented on the figure 
let s xy ^ 4/3, then 

(2 ~ S xy ) 


2A- 


A, 




and the stability area is given by (27) (figure 23). This area is reduced to V = 0 


when s xy = 2. 

Case 3 : s„„ < 2 s„. 


The inequalities (28) read 


{\±V x )(\±V y ) ^ 0, 
that is equivalent to (24). The inequalities ( |29[ ) write 

{\±V x )(\^±V y ) ^ 0, 


(31) 


for 7 = s xy /(2s q —s xy ). The inequalities (|3lj) have solutions if 7 ^ 0, that is equivalent 
to s xy ^ 0. If s xy ^ s , (31) contains (1241), that correspond to the same constraints 
that those imposed by (28). If -s xy ^ s q , then the stability area reduces to (26). 


There is still to study the influence of (30) on the stability area. These identities 
read 
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Figure 24. Layout of the areas given by (26) and (32). 


after division by 2 s g — s xy . We change of frame, V x = V x + V y , V y = V x — V y , and 


obtain 


4A 2 


2 ^ S X y 


(4 - (2s q + s xy )) ± 4 7 AV* + (V x ) 2 - (y y f > 0, 


and the analogous inequality exchanging the roles of V x and V y . We center these 


equations to obtain (25a) and (25b). 


Sub -case 1 : s xy ^ s q {2 — s q ). 


The inequality (25a) reads 


(V y ) 2 - ( V x ± 2A 7 ) 2 < 


16A 2 


(2Sq - S xy ) 2 


(s q (2 s q ) s xy ). 


(32) 


We represent the areas corresponding to the inequations (26) and (32) on the figure 
24 The relation (26) is equivalent in the new frame to |V| [ ^ 2A 7 . It is contained in 
the area delimited by (32) (figure 24). The final stability area is then given by (24) 
if 7 ^ 1 and by (26) otherwise. 


Sub -case 2 : s xy ^ s q {2 — s q ). 


The inequality (25a) has a solution V if the abscissa of the left pole of the hyperbole 

sgative 

2A 7 — 


centered in 2A 7 is nonnegative. This abscissa, equal to 

4A 

_I c _ 




(% y ^ 7 )) 2 1 


■xy) 


is nonnegative if 


12 s xy | ^ 2\s q 1|. 


(33) 


Assuming that s xy ^ 2, then s xy 


^ 2\s q — 1| + 2 ^ 2 s q , that gives an empty L°° 
stability area in V (refer to the beginning of the proof). As a consequence s xy ^ 2. 
If s q ^ 1, the inequality (33) is equivalent to s xy ^ 2 s q , otherwise to s xy ^ 2(2 — s q ). 
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Figure 25. Layout of the bounds of the L°° stability areas of the 
scheme relative to 2 = V with a non intrinsic diffusion for 7^1 and 

s 9 ( 2 - s q ) 


^ s xy ^ 1 - 


These conditions carrying on s ensure the existence of a non empty stability area in 

V. 


To obtain the final stability area in V given by the set of inequations from (28) to 
(30), we compare those associated with (25a) to (24) if 7 ^ 1 or (26) if 7 ^ 1. The 
two areas (24) and (26) read respectively ^ 2A, and 


|V |! ^ 2A 7 . 


(34) 


in the frame (0, V x ,V y ). If 7 ^ 1, the stability area is included in the one given by 
(34): it is represented by the points A, B , C , D on the figure 26 If 7 ^ 1 and s xy ^ 1 , 
tne area (24) is included in the one given by (25aJ: the abscissa of the left pole P 
of the hyperbole centered 2A7 is greater than 2A (figure 25). If s xy ^ 1 , it is not the 
case any more and the area decreases when s xy increases. □ 


The domain of L°° stability in s is identical for u = 0 and u = V. It is given by 
the triangle of the figures [TT] and [20] However, the choice u = V increases the set of 
the s verifying the optimal stability area \V loo < A: indeed the area ( |24[ ) spreads to 
TRT schemes consistent with s q ^ s xy ^ min(l,2s ? ). 


We now compare the areas in V associated with u = 0 and u = V for three s that 
don’t verify s q ^ s xy ^ min(l, 2s q ). We want to know if u = V provides a better 
L°° stability behaviour than 2 = 0. The associated inequalities are presented in 
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Figure 26. Layout of the bounds of the L°° stability areas of the 
scheme relative to u = V with a non intrinsic diffusion for 7^1 and 


Sq) ^ S X y 


Choice of s 

(0,1,1,1/2) 


(0,1,1,3/2) 




(0,3/2,3/2,3/4) 



Area for h=0 

(V x ±4) 2 —(V y ) 2 

7 12 

(V x ±4/3) 2 - 

-(V y ) 2 

> 

4/9 

(V x ±4) 2 

-(V y ) 2 


44/3 


(yy±4) 2 —(yz) 2 

> 12 

(V y ±4/3) 2 - 

-(V x ) 2 

;> 

4/9 

(V y ±4) 2 

-(V x ) 2 


44/3 

Area for u=V 

\V\eo < 1/3 


(V x ±6) 2 - 

(V y ) 2 

> 

32 

IH»< 

1/3 






(V y ± 6) 2 - 

(V x ) 2 

> 

32 








in»<i 









TABLE 3. L°° stability areas of the schemes with a non intrinsic dif¬ 
fusion for A = 1. 


the table [3] The solutions V are represented on the figures from 27 to 29 in the 
frame (O, V x , V y ): the areas delimited by the points A,B,C,D are associated with 
the scheme relative to u = V, and those bounded by E, F, G, H correspond to the 
MRT scheme (u = 0). 

We can not say that a scheme is better than the other in terms of L°° stability: 


generally the areas just intersect. In the section 2.3 we saw that the velocity held 


u = V improve the stability associated with the blow up of the scheme that rather 
corresponds to a L 2 notion than to a L°° notion. We confirm this fact theoretically 
in the section [4] 
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Figure 27. L°° stability areas for s = (0,1,1,1/2). Scheme relative 
to u = V (A, B,C, D), to u = 0 (E, F,G, H), with a non intrinsic 
diffusion. 



Figure 28. L°° stability areas for s = (0,1,1,3/2). Scheme relative 
to u = V (A, B,C, D), to u = 0 (E, F,G, H), with a non intrinsic 
diffusion. 
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Figure 29. L°° stability areas for s = (0,3/2,3/2,3/4). Scheme 
relative to u = V (A, B,C, D), to u = 0 (E, F,G, H), with a non 
intrinsic diffusion. 


3.2. The intrinsic diffusion case 

As for the non intrinsic case, we first deal with the MRT scheme corresponding to 

u = 0 . 

Proposition 3.3 (L°° stability areas for the MRT scheme). Let V £ M 2 , (s q . s xy ) £ 
M 2 , consider the twisted D 2 Q 4 MRT scheme 
parameters ( 0 , s q , s q , s xy ), and the equilibrium 

m eq (0) = p(l, V x , V y , 0). 


u = 0 ) associated with the relaxation 

m 


Note 7 = s xy /s q , 


• */ 0 ^ s xy < min (^, 2 - s q ) 
stable for all V such that 


(area BCD on the figure 30), the scheme is L c 


|V|i< a 7 . 


(35) 


if s q ^ min(l, s xy ) and s xy ^ 2 s q 
L°° stable for all V such that 


(area ABC on the figure 30), the scheme is 


V|i<A(2- 7 ). 


(36) 
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Figure 30. L°° stability area in s of the twisted D 0 Q 4 MRT scheme 
(u = 0) with an intrinsic diffusion. BCD : (35), ABC : (36), ACD : 

pit 


if s q ^ max(l, 2 — s xy ) and s xy ^ 2(2 — s q ) (area ACD on the figure 30), the 
scheme is L°° stable for all V such that 


|V| 1 <A(--2- 7 ). 


(37) 


if s xy = s q = 0 (point B on the figure 30), the scheme is unconditionally L c 
stable. 

For all other s, no V corresponds to a L°° stable scheme. 


In particular, the parameters ( s q ,s xy ) corresponding to a non empty area of L°° sta¬ 
bility in V are included in the square [0, 2] 2 . 


Proof. The positivity of the matrix sending f on f* reads 




2s n “f" SL-, S n 

q A y ) ± ^~(v x ± v y ) ^ 0 . 

4 4A 


(38) 


If s g = 0, then the inequalities (|38|) impose s xy = 0 and there is unconditional stability 


in V. Otherwise these inequalities are equivalent to 


|V|i ^ min (A 7 , A (2 — 7 ), A(- 2 - 7 )). 


(39) 
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The conditions to have a non empty set of V are to be determined: for that we must 
study when the minimum in (39) is nonnegative. 


7 ^ 0 <= 

o' 

A\ 

If 

7^07= 

^ S xy F 2 s y , 

7^07= 

^ S xy ^ 2(2 


We now suppose 
according to the 
s q < min(l, s xy ), 
closes the proof. 


that these conditions are verified and we determine the minimum 
choice of s: 


if 


^ minfs , 2 — s ), the area is given by (35), if 
it is given by (36), if s ^ max(l,2 — s x 


xy> 


we obtain (37). 


'his 

□ 


The areas are represented on the figure [30j We can deduce the L°° stability areas for 
the BGK scheme. 

Proposition 3.4 (The BGK case). Let V 6 M 2 , s£l, consider the twisted relative 
velocity D 2 Q 4 scheme, BGK of relaxation parameters s, and of equilibrium given by 

m 

• if 0 ^ s ^ 1, the scheme is L°° stable for all V such that 

|V|i < A. 

• if 1 ^ s ^ 4/3, the scheme is L°° stable for all V such that 

|V|i«A(t-3). 

• For all other s, no V corresponds to a L°° stable scheme. 

In particular, the parameter s corresponding to a non empty area of L°° stability in 
V are included in [0, 2]. 


We obtain a constant area of stability for s ^ 1. When s becomes larger than one 
(overrelaxation) the area decreases as s increases. For s greater than 4/3, there is no 
velocity L°° stable. Now we do the same job for the TRT scheme with u = V. 


Proposition 3.5 (L°° stability areas for a relative velocity scheme). Let V € 


' '%l 1 hi:y 


consider the twisted D 2 Q 4 scheme relative to u = V associated with 


the relaxation parameters ( 0 , s q , s q , s xy ), and the equilibrium (£| ~8) 


m eq (V) = p(l,0,U,-V x V y ). 
Note = (2 s q - s xy )/ ( s q - s xy ) and y 2 = s xy /{s q - s xy ), 
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• s q < s xy ^ min( 2 s ? , 2(2 
stable for all V such that 


s„)) (area ACD on the figure 31), the scheme is 


(V*±A 7l ) 2 -(V y ) 2 0 2 7l72 , 

(' V x ± A 71 ^ 72 ) 2 - (V y ± A) 2 < A 2 7 2, 

(y x ± - ( v y t a ) 2 < a 2 7 2 , 

(V x ± A 72 ) 2 - ( V y ) 2 < * 2 (4s 2 - 2 s q s xy - s 2 xy + 8 ( 5 ^ - s q )), 

\ S q S xy) 


(40a) 

(40b) 

(40c) 

(40d) 


plus the analogous inequalities exchanging V x and V v . 

if s xy ^ min(, 2(2 — s )) (area ABD on the figure 31), the scheme is stable 
for all V such that 


{V x ± A 7l ) 2 - (V y ) 2 ^ A 2 7l72 , (41a) 

C V x ± A 71 y 72 ) 2 - (V y ± A) 2 ^ A 2 72 , (41b) 

(V x ± A^^) 2 - (' V y 4= A) 2 ^ A 2 7 2 , (41c) 

(V x ± A 72 ) 2 - (V y ) 2 ^ , 2 (4s 2 - 2s q s xy - s 2 xy + 8 (s xy - s q )), (41d) 

\ s q s xy) 


plus the analogous inequalities exchanging V 
if s xy = s q (segment [AD] on the figure 31), the L 
the proposition \3.J\ 

For all other s, no V corresponds to a L°° stable scheme 


x and V y . 

stability area is given by 


In particular, the parameters ( s q ,s xy ) corresponding to a non empty area of L°° sta¬ 
bility in V are incliLded in the square [0,2] 2 . 


We now compare the L°° stability areas of the schemes relative to the velocities u = 0 
or u = V with an intrinsic diffusion for different s. We draw the areas given by the 


propositions 3.3 and 3.5 on the figures 32 and 33 


The results are similar to the ones obtained in the non intrinsic case. We can’t 
distinguish a scheme from an other in terms of L°° stability. The MRT scheme is 


better on the figure 32 but on the figure 33 the areas just intersect. 


4. Weighted L 2 stability 

In this section, we present some theoretical weighted L 2 stability results for the rela¬ 
tive velocity D 2 Q 4 schemes. These results, based on the notion of stability proposed 
by Yong and al in J3], confirm the phenomena numerically observed in the section 
In particular, taking u equal to the advection velocity V (“cascaded like” scheme) 


2.3 
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FIGURE 31. L°° stability area in s of the twisted D 2 Q 4 scheme relative 
to u = V with an intrinsic diffusion. ACD : (40a) to (40d|, ABD : 
dilal) to (Eldb. 



Figure 32. L°° stability areas for s = (0,1,1,3/2). Scheme relative 
to u = V (A, B, C , D ), u = 0 (E, F. G, H), with an intrinsic diffusion. 
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Figure 33. L°° stability areas for s = (0,3/2,3/2,3/4). Scheme 
relative to u = V (A, B, C, D, E, F,G, H), u = 0 (/, J, K, L), with 
an intrinsic diffusion. 


provides good stability features. These results also extend the L°° stability ones 
(section pi): the L°° stability areas are included in the L 2 ones. 


4.1. The weighted L 2 stability notion 


The general framework of the relative velocity D ( /Q q schemes, d, q E N*, is chosen to 
present the stability notion introduced in [3] and studied in mm- 


An iteration of a lattice Boltzmann scheme splits into a transport step plus a relax¬ 
ation step that is here linear since the equilibrium is linear ( ]6|7|8| ). Given the 

matrix composed of the distribution vectors taken on all the lattice points, we have 


f(.,t + At) = TR(u)f(.,t), t G M, 

where R(u ) is the relaxation matrix and T the linear non local transport operator. 
The size of the matrix is equal to the number of velocities q multiplied by 

the number of lattice points. In the case of the relative velocity D c iQ q scheme, the 
collision matrix reads R(u) = I q + J(u) where J(u) is given by 

J(u) = M(u)~ l DM(u){B - I q ), 

where D = diag(s) is the diagonal matrix of the relaxation parameters and B defines 
the linear equilibrium thanks to / eq = Bf. 


The calculus of the spectrum of the amplification operator TR(u) being difficult 
because of the matrix size, the idea is to study separately the transport and the 
collision: we want T and R(u) to be bounded by one in a well-chosen norm. The 
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amplification operator TR(u ) is then bounded by one and the scheme is stable for 
this norm. In the following, we say that an operator is stable in a given norm if the 
norm of this operator is bounded by one. 

The idea proposed in |3j consists in weighting the L 2 norm defined by 

q-i , 

ll/| 2 = y = ( y o> •' • > Vq-i) 

j =o 

so that the transport is an isometry and the collision operator is stable in the weighted 
norm. The classical L 2 norm keeps the transport as isometric but the collision norm 
is difficult to evaluate. Introducing a weight allows to overcome this difficulty. First, 
we define a norm on M 9 depending on an invertible matrix P E _Ad g (M), 

\y\ P = \Py\ 2 i y£^ q - 


We can then define a norm for a matrix g of the size q multiplied by the number of 
lattice points with 

\y\ p,c = ( \9( x )\p) 2 > 

x£C 

where g(x ) is the column of g associated with the node x. The necessity to define such 
a non local norm is due to the non locality of the transport. The matrix P is chosen 
such that f PP is diagonal. Thus the transport is an isometry for the norm | • | p£ 
for the periodic and bounceback boundary conditions |21j . Indeed, the hypothesis 
of “quasi orthogonality” fPP diagonal) carrying on P cancels all the cross terms in 
the calculation of |/(., t)\ p Then the isometry of the transport is obtained thanks 
to a simple change of variable: it uses the bijection between the nodes of the mesh 
and the nodes after the transport at a given velocity. 

Contrary to the transport, the collision is a local operator so that its study reduces 
to the vectorial norm | • | p . If the matrix R(u) is stable in a particular node x of the 
lattice C for | • | p , it is stable for | • | P £• That’s why, we use the local operator norm 

\\R(u)\\ p = sup \R(u)x \ p , 

\x\p=l 


for the collision. 


The matrix P is requested to diagonalize the collision: ||i?(u)|| p is then easy to 
evaluate. These requirements define a notion of stability structure first evocated in 
®- 

Definition 4.1 (|25|). A matrix N G _Ad g (M) has a pre-structure of stability if exists 
an invertible matrix P E At g (M) and some vectors s,p£l 5 so that 

PNP 1 = -diag(s), 
l PP = diag(p), 
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where diag(p) is the diagonal matrix whose diagonal is constituted of the coefficients 
of p. The pre-structure of stability becomes a structure of stability if 

e [0,2], 0O<g-l. (42) 


Suppose that J(u ) has a pre-structure of stability and consider the collision in norm 
| • |p. Knowing that the eigenvalues of R(u) are 1 — s k for 0 ^ k ^ q — 1, 

|P-R(u)P _ 1 Pcc | 2 = || PR( u)P 1 || 2 = max(|l — s|), 


I«(«)IIp = 


sup 
|.P:r | 2 =1 


where || • || 2 is the operator norm associated with | • | 2 . Thus the collision is stable for 
|| • ||p if the condition (42) is verified. Under these conditions, the scheme is stable 


in the norm 


I P,C 


since t 


re transport is isometric. 


We present a theorem from (25j giving a necessary and sufficient condition of existence 
of a pre-structure of stability. In the following, this theorem is the tool used to obtain 
some stability results for our relative velocity schemes. 


Theoreme 4.1 ([25]). A matrix N E _A/J g (M) has a pre-structure of stability if and 
only if there exists a diagonal positive definite matrix A E A4 g (M) so that 

NA = A*N. (43) 


This criteria gives a practical criteria of existence of a pre-structure of stability 
through the resolution of a linear system in the coefficients of A. The size of this 
linear system is q(q — l)/2 since the matrix IVA — A f IV is antisymmetric. 


The matrix P defining the weighted norm is explicitely derived in the proof of this 
theorem [25]. We exhibit P through the proof of the sufficient condition of pre- 
structure. The identity (43) is equivalent to the fact that A _1 IVA is symmetric with 
respect to the scalar product 

9-1 

(*> V)a = A u*> y e Rq - 
* ,1=0 


This implies the existence of an orthonormal matrix Q in the sense of (-,-)a diag¬ 
onalizing A - 1 IVA. Then, the matrix P = (AQ )” 1 diagonalizes the collision. The 
matrix t PP is diagonal because of the orthonormality of Q for the scalar product 


4.2. Stability results for the twisted relative velocity D 2 Q 4 scheme 


In this section, we apply the stability criteria given by the theorem 4.1 to obtain 
weighted L 2 results for the twisted relative velocity D 0 Q 4 scheme. Our motivation is 
to compare the cases u = 0 and u = V. To do so, the scheme must have at least 
two different relaxation parameters because the BGK scheme does not depend on the 
relative velocity since ([3]) is verified. We begin by studying the MRT scheme (u = 0) 
whose results are limited. We then obtain more results for the cascaded like scheme 
(u = V). 
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4.2.1. The MRT scheme 

We focus on the MRT scheme corresponding to u = 0. We limit to V = 0 because 
when u = 0 and 7/0, there is no pre-structure of stability. The framework does 
not distinguish the non intrinsic and intrinsic case because for V = 0, the equilibria 
Q and ([8]) are identical. 

Proposition 4.1 (L 2 stability for the MRT scheme). Consider the twisted D2Q4 
scheme relative to u = 0 , of equilibrium m eq (0) = (p, 0 , 0 , 0 ) associated with the 
relaxation parameters s = (0, s q , s q , s xy ) E M 4 . The associated matrix J( 0) has a 
pre-structure of stability. Moreover, if 0 ^ s q , s xy ^ 2, then J( 0) has a structure of 
stability. The collision matrix R( 0) is then stable for || • ||p and the scheme is stable 
in norm I • I p r . 


Proof. Since u = V = 0, the matrix J(0) reads 

J(0) = M(0)~ 1 D(E - J 4 )M(0), 

with the diagonal matrix E = diag(l, 0,0,0). The matrix D(E — J 4 ) is then diagonal 
as a product of two diagonal matrices. Since M{ 0) -1 = f lW(0)/4, the matrix 

J(0) = t M(0)D(E — / 4 )AT(0)/4 is symmetric. So the identity is solution of the 
equation (43) and «/(0) has a pre-structure of stability. The spectrum of J(0) being 
composed by 0, — s g , — s g , — s xy , the result on the structure of stability is obvious. □ 


4.2.2. The scheme relative to the advection velocity (“cascaded like” scheme) 

We now focus on the choice u = V beginning by the non intrinsic case. 

Proposition 4.2 ( L 2 stability for a non intrinsic relative velocity scheme). Consider 
the twisted D0Q4 scheme relative to u = V = (V x ,V y ) E M 2 , of equilibrium 

m eq (V) = p(l, 0,0,0), 

associated with the relaxation parameters s = (0, s q , s q , s xy ) E M 4 . The matrix J(V) 
has a pre-striLctiLre of stability if and only if 1^1^ < A. Moreover, if 0 ^ s q ,s xy ^ 2, 
then J(V) has a structure of stability. The collision matrix R(V) is then stable for 
|| • ||p and the scheme is stable in norm \ ■ |p £ - 


Proof. According to the theorem 4.1, the existence of a pre-structure of stability for 
N = J(V) is equivalent to the existence of a diagonal positive definite matrix A 
such that 

J(V )A = A f J(F). 

A solution of this linear system is given by 
f (X + V x )(X + V y ) 0 

(X~V X )(X + V y ) 


A = 


V 


0 

0 

{X-V x ){X-V y ) 

0 


(.X + V x ){X-V y )J 
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This matrix is positive definite if and only if IV^ < A, that closes the first part 
of the proof. The eigenvalues of J(V) are 0 , —s q ,—s q ,—s xy because the matrix 
M(0)J(V)M(0)~ 1 is given by 

/ 0 0 0 0 \ 
s q V x —s q 0 0 

SqV y 0 — Sq 0 

\(2 Sq -S xy )V*Vy vy(s xy -s q ) V X (s xy -Sq) sj 
We deduce that 

PJ(V)p- 1 = diag(0, -s q , -Sq, -S xy ), 


and that J(V) has a structure of stability if 0 ^ s q , s xy ^ 2. The norm || • ||p of the 
collision operator R(V) is equal to 1 and the scheme is stable in norm | • |p □ 


proposition 
weighted P 


3.2 


These results extend the L°° stability ones obtained in the proposition |3.2| The 
gives the L°° stability area | ^ A for s q ^ s xy ^ min(l, 2s ). The 

u, xy SC 2. This result was expected because 


notion generalizes it to 0 ^ s x 


the numerical experiments of the section 2.3 put in evidence the independence of the 
area of L 2 stability towards the relaxation parameters. The fact that we can not 
obtain the same type of result for the MRT scheme as for u = V is an other evidence 
of the good behaviour of this latter. 


Let’s also note that these stability conditions are defined by opened sets. The nu¬ 
merical experiments seem to show that the scheme is still L 2 stable on the closure 
of these sets. However, it is not possible to proove it with this notion because the 
matrix A is null on this closure. This proposition leads to a natural corollary for 
a single relaxation parameter (BGK) that has already been proved in |25]. In the 
BGK case, the scheme does not depend on u. In particular, MRT and cascaded like 
scheme are identical. That’s why the following result is valid whatever the relative 
velocity u. 

Corollary 4.1 (The BGK non intrinsic case [25]). For V = ( V x , V y ) E R 2 , consider 
the twisted relative velocity D2Q4 BGK scheme of equilibrium 

m eq (u) = p( 1, V x - u x , V y - u y , {V x - d x )(V y - u y )), 

associated with the relaxation parameter s £ M. The matrix J(u ) has a pre-structiLre 
of stability if and only if | < A. Moreover, if 0 ^ s ^ 2, then J(u ) has a structiLre 

of stability. The collision matrix R(u ) is then stable for || • || p and the scheme is 
stable in norm \ ■ |p £ - 

Proposition 4.3 (L 2 stability for an intrinsic relative velocity scheme). Consider 
the twisted D2Q4 scheme relative to u = V = (V x ,V y ) £ R 2 , of equilibrium 

m eq (V) = p(l, 0,0,-V x V y ), 

associated with the relaxation parameters s = (0, s q , s q , s xy ) £ R 4 . Suppose that V x = 
0 (resp. V y = 0) then the associated matrix J(V) has a pre-structure of stability if 
and only if \V V \ < A (resp. \V X \ < X). Moreover, if 0 ^ s q ,s xy ^ 2, then J{V) has 
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a structure of stability. The collision matrix R{V) is then stable for || ■ || p and the 
scheme is stable in norm I • I p r . 


Proof. A non null solution of the equation (43) exists if and only if one of the two 
components of V is null. Suppose that it is V y . then a solution of (43) is given by 

/\ + V x 


A = 


V 


o 

o 

0 


0 

X-V x 

0 

0 


0 

0 

A-K* 

0 


This matrix is definite nonnegative if and only if \V' 


0 

0 

0 

A + K* 
c | < A. 


structure of stability is identical to the one of the proposition 4.1 


\ 


The reasoning to get a 

□ 


This proposition extends the L°° stability results of the proposition 3.5 for the di¬ 
rections V x = 0 and V y = 0. The L°° notion provides areas decreasing as much 
as s q and s xy go far from each other. The | • | p/ - notion gives a constant optimal 
area in V for these directions while these parameters are bounded by 0 and 2. This 


phenomenon was observed numerically on the figures 11 to 16 


We close the section by a proposition for the BGK case that is not a corollary of the 
proposition for the intrinsic cascaded like scheme. 


Proposition 4.4 (The BGK intrinsic case). For V = ( V x ,V y ) E M 2 , consider the 
twisted D 2 Q 4 scheme, of equilibrium 

m eq (S) = p( 1, V x - u x , V y - u y , u x u y - u x V y - u y V x ), 

BGK of relaxation parameter s E M. The associated matrix J{u) has a pre-structure 
of stability if and only if\V\ l < A. Moreover, ifO ^ s ^ 2, then J(u ) has a structure 
of stability. The collision matrix R(u) is then stable for || • || p and the scheme is 
stable in norm I • I p r . 


We note tha t th is proposition generalizes the L°° stability results in the BGK case 
(proposition |3.4[ ) . This weighted L 2 notion extends the set of s corresponding to the 
stability area < A to 0 ^ s ^ 2 . 


4.2.3. Interpretation of the results 


When we compare these propositions to the numerical results obtained in the section 
2.3, we notice that this stability notion is only usable when the parameters s and V 


are decoupled. For the scheme relative to u = V with an intrinsic diffusion, we have 
only theoretical results for V x = 0 ou V y = 0, corresponding to a case where the 
areas in V are independent of s. Instead, when the area in V is a function of s, it is 
not possible to build a pre-structure of stability any more. For example, there is no 
pre-structure of stability when V x and V y are different from 0 for the D 2 Q 4 scheme 
relative to u = V with an intrinsic diffusion: for this scheme, the numerical results of 


the section 2.3 exhibit a link between V and s. Similarly, the MRT scheme present 
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stability areas in V depending on s and it is impossible to build a pre-structure of 
stability for V / 0. 


The origin of this limitation seems to be the hypothesis of diagonalization of the 
collision. Indeed, the spectrum of this operator being equal to 0, — s , — s , — s xy , the 
existence of a pre-structure of stability with such a diagonalization implies automat¬ 
ically the stability of the scheme for 0 ^ s , s xy ^ 2. Then V and s can’t be linked 
because the former hypothesis is too requiring. The linear system of six equations 
and four unknowns for the D 2 Q 4 scheme is an evidence of the constraints imposed 
by the notion. The existence of a pre-structure of stability requires its rank to be 
three that imposes V and s to be decoupled here. For example, when u = 0 and 
F/ 0 , there is no solution positive definite to the equation |43[ ), because this rank 
is greater than four. That’s why it is not possible to show similar results than in the 
case u = V which requires the resolution of a rank three linear system. 


However, this notion gives some promising results as those exposed here and in 
eh Eg. It seems to describe well the limit of blow-up of a scheme. The main purpose 
now is to obtain theoretical stability results in more complex cases. Particularly, we 
want to complete this study for the D 2 Q 4 especially when u = 0. To do so, we must 
be able to relax the constraint of diagonalization of the collision without penalizing 
the isometry of the transport. Maybe, it would be possible to find some matrix P so 
that the norm • Ip £ of the collision operator is computable without diagonalizing 
it. 


Conclusion 

We have studied the stability of a four velocities relative velocity scheme with two 
different equilibria for a linear advection equation. The discussion is based on two no¬ 
tions of L°° and weighted L 2 stability. It carries on the choice of the relative velocity 
equal to 0 (MRT scheme) or to the advection velocity (“cascaded like” scheme). The 
main conclusions are the following: comparing MRT and relative velocity schemes, 
no scheme is “better than the other” in terms of L°° stability; the stability areas 
generally just intersect. Instead, in terms of L 2 norm, the scheme relative to the 
advection velocity is more stable than the MRT scheme. This improvement is cor¬ 
related to the cancellation of some dispersive terms of the equivalent equations for 
the scheme relative to the advection velocity. These terms stay for u = 0 and are 
linked to instabilities at low diffusion. This result has been justified theoretically for 
the scheme relative to the advection velocity thanks to the weighted L 2 notion of 
stability. Further work needs to be done to obtain more precise theoretical results 
for the MRT scheme and for systems of conservation laws. The study also needs to 
be generalized to an arbitrary relative velocity u, the present work being adapted to 
the two choices u = 0 and u = V. 
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Appendix A. Link between the twisted D 2 Q 4 and the D 2 Q 4 schemes 

We proove that the L°° and L 2 stability areas of the twisted D 9 Q 4 and D 2 Q 4 relative 
velocity schemes with a non intrinsic diffusion correspond through the composition 
of a rotation and a homothety. The same result is still true in the intrinsic case. 

If we note v the velocity set of the D 2 Q i scheme, the twisted set of velocities is given 
by 7 Zv where 

1 - 1 ' 


n = 


1 1 


is the composition of a rotation and a homothety of scale factor \/2. This transfor¬ 
mation allows to link the relaxation of both schemes. 

Lemma A.l (Relation between the relaxation operators). Let V £ R 2 , u £ R 2 , 
( s q ,s xy ) £ R 2 , the relative velocity D 2 Q 4 scheme associated with the equilibrium 

m eq (u) = p( 1, V x - u x , V y - u y , ( V x - u x ) 2 - ( V y - d y ) 2 ) , (44) 

for the relaxation parameters ( 0 , s , s , s ) and the relative velocity twisted D 2 Q 4 
scheme associated with the equilibrium 

m eq (u)=p(i, ( nv) x -{nu ) x , (• iiv) y -(nu ) y , 

((7ZV) x -(7Zu) x )((7ZV) y -(7Zu) y )), 

and with the relaxation parameters ( 0 , s q , s q , s xy ) have the same relaxation phase. 

Proof. We have defined a transformation 7 Z sending the velocities of the D 2 Q 4 scheme 
on the velocities of the twisted D 2 Q 4 scheme. This naturally leads to the following 
application also denoted by 7 Z. 


K : R[X, Y] -> R[A, Y] 

P eA 7^(P) ’ 


where 


7 Z(P) : 


M 


24 V 


V = { Vj , 0 ^ j ^ 3} HA P(fR,v j )f j • 

i=o 

Indeed the images of the moments 1, X, Y, X 2 — Y 2 by 7 Z are 

7*(1)(V) = 1(V), 

3 3 

7^(A)(V) = X(7^t;.)/. = 2(«/ - i)fj = ( X - y )0^ 

7=o 7=o 

3 3 

t^OO(v) = E Yinv^fj = + $)fj = ( x + Y )(V), 

7=0 7=0 


(45) 
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nxY)(v) = Y,( nv j) x ( nv j) y fj = E ((«/) 2 - (f) 2 )/j = (^ 2 - y2 )( v )- ( 46 ) 

3=0 3=0 

Providing these images, we can write the transformation of the relaxation of the D 2 Q 4 
relative velocity scheme by 1Z. We here choose to keep the polynomial notations to 
express the moments. The relaxation of the D 0 Q 4 reads 


X*(V) = X(V) + s q ((V x - u x )p - X(V)), 

Y*(V) = Y(V) + s q {(yy - u y ) P - y(v)), 

(. X 2 -Y 2 )*{V ) = {X 2 -Y 2 )(V) +s xy (((V x -u x f - (V y -u y ) 2 )p-(X 2 -Y 2 )(V)]. 


Introducing the relations (45|46), we obtain 


iz(xnv) = nx)(v) + s q {((KV) x - (Ku) x )p - n(x)(v)), 

K(Y)*(V) = n(Y)(V) + s q {{(KV) y - (Ku) y )p - U(Y)(V)), 

K(XY)*(V) = K(XY)(V) + s xy (((nV) x -(Ku) x )«KV) y -(Ku) y )p-K(XY)(V)), 


that is the relaxation of the twisted D 2 Q 4 relative velocity scheme. Note that this 
last step uses the fact that X and Y have the same relaxation parameter s q . □ 


Proposition A.l (Relation between the stability areas). Let (s q , s xy ) £ M , note 
$nt C M 2 (resp. S t C the set of the velocities V £ M 2 such that the relative 
velocity D 2 Q 4 scheme (resp. twisted relative velocity D 2 Q 4 scheme) of equilibrium 


given by (44) (resp. and of relaxation parameters s = (0, s q , s q , s xy ) is L' 


or 


L stable. We have 


S t — LlS nt . 


(47) 


Proof. Note R t (V,u, s) £ A4 4 (M), the relaxation operator of the twisted relative 
velocity D 2 Q 4 scheme of equilibrium (6|7): it verifies 

f* = R t (V,u,s)f. 

Note R nt (V,u, s) £ A4 4 (M) its analogue for the D 2 Q 4 scheme associated with the 
equilibrium ( |44[ ). The lemma |A.l means that 

Rnt( v > s ) = R t (KV, Hu, s ). (48) 


The transport does not influence the L°° stability because it only exchanges the 


V £ 


so that the 


particle distributions. Thus '-’rat is the set of the velocities 
matrix R nt (V,u, s) is nonnegative. According to the relation (48), it corresponds 
to the velocities V £ M 2 so that 7 ZV belongs to the L°° stability area of the twisted 


scheme: that is equivalent to the relation (47) and the lemma is proven in the L 


case. 
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Instead, the transport plays a role for the L 2 stability. It becomes local in the Fourier 
space: the transport matrices associated with the D2Q4 scheme (A nt ) and the twisted 
one (A t ) are given by 

A nt (k ) = diag(exp(ztr ■ k), 0 < j < 3), fc E M 2 , 

A t (k ) = diag(exp {iTZv- ■ k),0 ^ j ^ 3), k E M 2 . 

The transport matrix A nt {k) is equal to A t (7Zk/2) because TZv- ■ TZk/2 = v- ■ k 

since TZ is the composition of a rotation and a homothety of scale factor \[2. The 
amplification matrices L nt and L t of both schemes are then defined by 

L t (k, V, u , s) = A t (k ) R t {V, u, s), 

L nt {k, V , u , s) = A t {lZk/ 2) R t {TZV, TZu , s). 

We then have the following identity 

L nt (k, V , u, s ) = L t (lZk/ 2,7£V, TZu , s). (49) 


The L 2 stability needs to evaluate the maximum of the spectral radius r of this 
matrix for all the wavenumbers k E M 2 . According to (49), 

max r(L t (k,V,u, s )) = max r(L t (TZk/2, TZV, TZu, s)) 
km 2 fee R 2 

= max r(LJk, TZV, TZu , s)), 
fceR 2 

the last equality coming from a variable change in M 2 . As for the case of the L°° 
stability, only the velocities V finally matter: the relation between the L 2 stability 
areas is then obtained in the same way. □ 


Appendix B. Theoretical L°° stability results for the D 2 Q 4 scheme 

B.l. For a non intrinsic diffusion 

Proposition B.l (L°° stability areas for the MRT scheme). Let V € M 2 , (s q . s xy ) G 
M 2 , consider the D 2 Q 4 of MRT scheme with the relaxation parameters (0, s q , s q , s xy ), 
associated with the equilibrium 

m eq (0) = p{ 1, V\ {V x f - ( V v f ). 

Note 7 = s q /s xy , 

• if 0 < s xy ^ min(s g , 2 — s q ), the scheme is L°° stable for all V so that 

(' V x ± A7) 2 - (V y ) 2 > A 2 ( 7 2 - 1 ), 

( V y ± A7) 2 - (V x ) 2 ^ A 2 (7 2 - 1 ). 

• if s q ^ s xy ^ 2 s q and s q ^ 1, the scheme is L°° stable for all V so that 

(V x ± A7) 2 — (V y ) 2 ^ A 2 ( 7 — l) 2 , 

(■ V y ± A7) 2 - (V x ) 2 ^ A 2 (7 - l) 2 . 
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• if 2 — s q ^ s xy ^ 2(2 — s ) and s q ^ 1, i/ie scheme is L°° stable for all V so 
that 

{V x ± A 7 ) 2 - (V y ) 2 > A 2 ((7 + l) 2 - , 

xy 

{ v y ± a 7 ) 2 - ( f *) 2 > a 2 (( 7 +1 ) 2 - . 

xy 

• if s xy = 0 and 0 < s g ^ 2, f/ie scheme is L°° stable for V = 0. 

• if s xy = s q = 0, i/ie scheme is unconditionally L°° stable. 

• For all other s, there is no V corresponding to a L°° stable scheme. 

Proposition B.2. Let V £ M 2 , (s q , s xy ) £ M 2 , consider the D 2 Q 4 scheme relative to 
u = V associated with the relaxation parameters (0, s q , s q , s xy ), and the equilibrium 

m eq (V) = p(l, 0,0,0). 

Let us note 7 = s xy / (2 s q — s xy ), 


• if s q ^ s xy ^ min(l, 2 .^), the scheme is L°° stable on 

|V| X <A. (50) 

• if s xy < 2 s q , max(s ? , 1) ^ s xy ^ 2(2 — s ), the scheme is L°° stable on the 
intersection of with 

a \2 

(V* ± A 7 ) 2 - (V y ) 2 > _ (^ - ,,(2 - q,)), (51a) 

4A 2 

(V s ± a 7 ) 2 - cn 2 > , 9q _ g - v 2 -*))■ ( 51b ) 

• if 0 ^ s xy ^ min(sq,s ? (2 — s )), f/ie scheme is L°° stable on 


|V|!< A 7 . 

if s q (2 — s q ) ^ ^ min(s g , 2(2 — s ? )), i/ie scheme is L°° stable on ( 51a^51b ). 

if s xy = 2 s q and 1 < s xy ^ 2, f/ie scheme is L°° stable on the intersection of 


(50) and 


l^loo^A 


2 - a 


xy 


xy 


• */ s xy = = 0, i/ie scheme is unconditionally L°° stable. 

• For all other s, no V corresponds to a L°° stable scheme. 


B.2. For an intrinsic diffusion 

Proposition B.3. Let V £ M 2 , (s g ,s xy ) £ M 2 , consider the D 2 Q 4 MRT scheme 
associated with the relaxation parameters (0, s , s , s ), and the equilibrium 

m eq (0) = p (l,p a: j p2/,0). 

IVote 7 = 
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if 0 ^ s xy ^ min^, 2 — s q ), the scheme is L°° stable for all V such that 


l^loo^ 


A 7 


if s q ^ min(l, s xy ) and s xy ^ 2 s Q the scheme is L°° stable for all V such that 


xy 


IVL, £ 


A(2- 7 ) 


if s q ^ max(l, 2 — s xy ) and s xy ^ 2(2 — s q ) the scheme is L°° stable for all V 
such that 


IVL, £ 


A(i- 2 - 7 ) 


• s xy = s q = 0, the scheme is unconditionally L°° stable. 

• For all other s, no V corresponds to a L°° stable scheme. 

Proposition B.4. Let V € M 2 , (s q , s xy ) 6 M 2 , consider the D 2 Q 4 scheme relative to 
u=V associated with the relaxation parameters (0 ,s q ,s q ,s xy ), and the equilibrium 

m eq (V) = p ( 1, 0,0, {V y f - (V x ) 2 ). 

Note 7 j = (2 s q - s xy )/ (s g - s xy ) and ~/ 2 = s xy /(s q - s xy ), 

• if s q < s xy ^ min(2s g ,2(2 — s )), the scheme is stable for all V such that 

^7lx2 firy\2 ^ ^hl2 


(' V X db 771 )2 _ (yy)2 ^ 


(V x ± A ^^) 2 - (V»=F ^) 2 ^ 


( V x db A 


7! +72x2 


A 2 7 2 2 


4 - (V* ± -) 2 < 4 , 


(4 S q 2 S q S xy S xy + 8(s xy 


plus the analogous inequalities exchanging V x and V y . 

if s xy ^ min(s ? , 2(2 — s q )), the scheme is stable for all V such that 

(v* ± ^) 2 - (yyf ^ A ^ 172 , 


,7! +72x2 


(V x ±\ 


(v x =b a 11 ' A 12 ) 2 — (y y ip ^) 2 ^ " 72 


7i+7 2 n2 nrv x^2 


A 2 7 2 2 


4 >'-(^± 3 )^ 4 


4 ’ 

A 2 7 2 




(4Sg 2 S q S xy S xy + 8(s x y 


plus the analogous inequalities exchanging V x and V y . 
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if s xy = s q ^ 1, the scheme is L°° stable for all V such that 




if 1 ^ s xy = s q ^ 4/3, t/ie scheme is L°° stable for all V such that 


|VL^-3). 

For all other s, no V corresponds to a L°° stable scheme. 


Appendix C. Proof of the proposition 13.51 


The positivity of the matrix shifting from f to f* is equivalent to 

A(2s 9 - s xy )(\ ± ( V x + (-1 TV*)) + 2(-l Y(s q - s xy )V x V y > 0, (52a) 

A 2 ^ ± A (s xy V x + (~iy(2s q - s xy )V y ) + 2(-l y\s q - s xy )V x V y > 0, (52b) 

A 2 s xy ± A (s xy V y + (-lY(2s q - s xy )V x ) + 2(-l )\s q - s xy )V x V y > 0, (52c) 

4A 2 - {{2s q + s xy ) A 2 ± s xy X(V x + (-1 yv y ) - 2(-l )\s q - s xy )V x V y ) > 0. (52d) 


for i = 0 , 1 . 
proposition 


We begin by the case s a = s xy of a BGK scheme corresponding to the 


3.4 The four inequations (52a) to (52d) become 


Xs q (X ± (V x + (-l) l V y )) ^ 0, 

(53) 

Sq (x±(v x + (-iyv y ))>o, 

(54) 

Sq (x±(v y + (-iyv x ))>o, 

(55) 

- s g (3A 2 ± X{V X + (-1 yv y ) ^ 0 . 

(56) 


The case s q ^ 0 is impossible because no velocity V would verify these four inequali¬ 
ties requiring V x ± V y to be simultaneously greater than A and inferior to —A. Then 
s q ^ 0 that leads to 


\V\ 1 ^ min(A, A(-3)), 


that gives exactly the proposition 3.4 according to the choice of s q . 
can write the inequalities (52a) to (52d) thanks to and 7 2 . If s xy 


If s xy y s q , we 
< s q , we obtain 


\^{X±(V X + (-iyv y )) + (~iyv x v y >0, * = 0,1, (57a) 

x 2 ^±^ 2 v x + (-iy^v y ) + {-iyv x v y ^o, » = o, 1, (57b) 

A 2 ^±^( 7 2 ^ + (-1)Si^) + (-1)^P^0, i = 0,1, (57c) 

4 - ( 2fi g + s x y ) A 2 ± y*x{ V x + (-1 yv y ) + (-1 yv x v y ^ 0 , * = 0,1. (57d) 

2 (^ - s xy ) 2 
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We do the following change of variable V x = V x + V y , V y = V x — V y . Noticing that 
— 7 2 = 2, the inequalities (57a) to (57d) are equivalent to 

2A 7l (A ± V x ) + (V x ) 2 - (V y f ^ 0, 

2A 2 7 2 ± 2A( 7l + 72 v j - V y ) + ( V x ) 2 - ( V y ) 2 ^ 0, 


2A 2 7 2 ± 2A( 7l + 7a v x + V y ) + {V x f - ( V y ) 2 ^ 0, 
8 - (As q + 2s xy )_ a2 ± ^ xvx + (ya;)2 _ {V y )2 ^ 


Sq $xy 


(58a) 

(58b) 

(58c) 

(58d) 


for 1 = 0 in addition to the ones obtained exchanging V x and V y corresponding to 


1 = 1. We center these inequalities to get the identities from (41a) to (41d). The 


case s ^ s is obtained by changing the sense of the inequalities. It is characterized 


by the inequalities from (40a) to (40d). 


We begin by treating the case s q < s xy corresponding to the inequalities from (40a) to 
(40d). The case s q > s xy will be considered further. The reasoning first eliminates the 
couples (s ,s x ) leading to no stable velocity and then concentrates on stable areas. 


First, we assume that s > max(0, 2s ) and show that there is no stable velocity V. 


Since s > s , we have r ) 1 ^ 2 ^ 0 and the equation (40a) reads 


(v v ) 2 — (y x ± A 7l ) 2 ^ —A 2 7 l 7 2 , 

plus the two ones obtained exchanging V x and V v . The intersection of these four 
hyperbolic areas is empty. 

Secondly we eliminate the case s q < s xy < 0. In this situation, 7 1 7 2 is nonnegative 
( 7 X and 7 2 are both nonnegative) and (40a) has a solution if and only if the abscissa 
of the right pole of the hyperbole centered in —Aq x is nonnegative. This reads 


-A71 + A( 7 j 7 2 ) 2 ^ 0 


7i < 7 2 


$xy 7 S q . 


The last inequality contradicts our hypothesis and there is no stable velocity when 


7 S xy < 0. 


It remains to treat the case s q < s xy ^ 2 s q for which 7 1 7 2 ^ 0. The equation (40a) 
has some solutions V if and only if the abscissa of the right pole P of the hyperbole 
centered in Ay x is nonnegative 


A71 + A(7i7 2 )2 ^ 0. 


(59) 


The area satisfying these inequations corresponds to the points A, B, C, D on the 
figure [35] It increases as the abscissa of P increases. The inequality (59) is equivalent 


to 7 X ^ 0 , that is verified since s q < s xy ^ 2 s . 


For the equation (40d), there are two cases. Hence 


4.s„ 2 s q s xy s xy + 8(s xy s q ) ^ 0, 
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Figure 34. Hyperbole As q —2s q s xy — s xy +8(s xy — s q ) = 0, and straight 
lines s xy = 2 - s q and s xy = s q . 



FIGURE 35. Layout of the bounds of the L°° stability areas corre¬ 
sponding to (40a) and (40d) for the twisted D 2 Q 4 scheme relative to 
u = V with an intrinsic diffusion. 


then the reasoning is the same as the one made in the case s xy > max( 0 , s q ) to show 
that there is no stable velocity. Hence 

4-Sq — 2 S q S X y — s xy + 8 (s xy Sq ) ^ 0, 
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then the abscissa of the right pole of the hyperbole centered in Ay 2 must be nonneg¬ 
ative: 

^72 + — I ^SqS X y S X y + 8 (s x y 'Sq))^ ^ 0. 


I $ X y | 


The parameter y 2 being negative, this is equivalent to 


1 


(4s 2s s s x y T 8(s ^q)) ^ 72> 


and so to 


since s q ^ s xy . 


( s q ~ s xyY 

(s q - ^ y )(4 s q + 2 s xy - 8) ^ 0 <F=> s xy s: 2(2 - s g ), 


We now summarize the results we have obtained for s q < s xy . We have proven that if 
s q < s xy . the inequalities (40a) and (40d) have solutions for if s xy < min(2s g ,2(2-s q )). 
It remains to prove that (40c) and ( |40b ) have also solutions in this case. 

The conditions (40c) and (40b) lead to a non empty area in V since the intersection 
I (figure 36) between the hyperbole given by (40c) centered in (—A(q x + 7 2 )/2, A) 
and the straight line V y = 0 has a negative abscissa. The same condition on the 
intersection of the hyperbole centered in (—A(7 1 +7 2 )/2, —A) and V y = 0 is requested 
for the equation (40b). We focus on the hyperbole centered in (—A(y x + 7 2 )/2, A) 
associated with (40c), because the one associated with ( |40b ) leads to the same result. 
The intersections of the hyperbole of equation 


(V x + A 


7! +72x2 


) z — (y v — x) z = a 2 7 2 , 


with V y = 0 verify 


(V- + A^p2) 2 - A 2 = A 2 


2 

72- 


Their abscissas are then given by 


-a^±a (7 = + i)L 


We check when the abscissa of the point I of the figure 36 is negative 

'h + 7 2 \ 2 


_ A 7p+72 _A( 7 2 + i)^o 


<72 + 1 


2/2 
s q < s xy 


+ ( s q ~ s xy ) 


This is equivalent to s xy {s xy — s q ) ^ 0, that is valid because s xy ^ s q . 

We now focus on the case s xy < s q corresponding to the inequalities from (41a) to 
(41d). Let’s begin by eliminating the case 2 s q ^ s xy < s q . In this case (41a) has 
solutions if and only if both poles of the hyperbole centered in —Ay x (y x ^ 0) have 
nonnegative abscissas that reads 


-A(7i + (7i7 2 ) 2 ) > 0- 


Because of the negativity of 7 1 , this is equivalent to s q ^ s xy that contradicts our 
hypothesis. 
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FIGURE 36. Layout of the bounds of the L°° stability areas corre¬ 
sponding to (40c) and (40b) for the twisted D 2 Q 4 scheme relative to 
u = V with an intrinsic diffusion. 


Second, we eliminate the area 2(2 — s ) ^ s xy < s q . This is done showing that there 
is no solution to (41d). In this area, we have 


4 — 2 <? 

^Oq A&q& X y 


$xy "F 5^) ^ 0, 


■xy 


(figure 34) so that (41d) has a solution if and only if the abscissas of both poles of 
the hyperbole centered in A|q 2 | are nonnegative. This reads 

A 


a It 2 I - 


Sq S X y | 


( 4 Sq ~ 2s q S xy ~ S xy + 8(s_ - %))5 ^ 0 


that is equivalent to s xy ^ 2(2 — s q ). This contradicts our hypothesis. 

It remains to eliminate the area s < min(0, s q , 2(2 — s )) corresponding to ^ 0 


and 7 2 ^ 0. If s q ^ 0, then 


7i +7 2 


0 and the inequality (41c) is non empty if the 


1 " ~ 7 - 2 

abscissas of the intersections of the hyperbole centered in (—A(y x + 7 2 )/2, A) with 
V y = 0 are both nonnegative. This is equivalent to 


a 7l±72 


A( 72 2 + 1)^0 


7i +7 2 


> 72 + 1 


^ S ly + ( S q 



This is equivalent to s xy (s xy — s q ) ^ 0, that is verified for s xy ^ s q since s xy < 0. This 
enters in contradiction with our hypothesis. If s > 0, the reasoning is symmetric 
with respect to the ordinates axis ( V x = 0) and leads to the same contradiction. 
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To close the proof, it remains to guarantee the existence of at least one velocity stable 
in the triangle defined by 0 ^ s xy ^ min(Sq, 2(2 — s )). The equation (41a) has a non 
empty stability area if and only if the two poles of the hyperbole centered in A 7l have 
a nonnegative abscissa: this case, illustrated by the figure [T9l is equivalent to 


Kh ~ (7i7 2 ) 2 ) > 0. 

This inequality is equivalent to 7l ^ 0, that is true for s xy ^ s . 


For the equation (41d), we distinguish two cases: the first is 


4 si - 2 - s 2 + 8(s - s ) ^ 0. 


9 9 xy "it/ ' “\“xy q) 

There is a non empty stability area if and only if the two poles of the hyperbole 
centered in A 7s have nonnegative abscissas (represented on the figure [l9|) : this reads 

^7 2 - i \ i ( 4 ^ - 2 s Q s xy - s 2 xy + 8 (s xy - q,)) 2 ^ 0. 


I Sq S xy | 


This inequality is equivalent to s xy ^ 2(2 — s q ). Otherwise, 

" 2 " - 4y + 8(s~, - s a) < 0, 


4 s — 2s ‘I 


and the equation (41d) reads 

(V y ) 2 - (V x ± a 72 ) 2 ^ - 


A 2 


( $xy ) 


2 (4Sg 2SgS xy Sjy + 8(Sj,j ^g))- 


This inequation has always solutions: the area of stability in V is analogous to the 
one represented on the figure [24] 

The inequations ( 41c[ ) and ( |41b[ ), illustrated on the figure [37j have solutions since the 
two intersections of the hyperbole centered in (A( 7l + 72 )/2, —A) with V y = 0 have 
nonnegative abscissas 


a 2i+ 72 _ A(7 2 + jjl o 


Tl + TA% 7 | + 1 


2 \ 2 
Sq ^ S xy 


+ (S„ - S xy ) 2 . 


This is equivalent to s xy (s xy — s q ) ^ 0, that is verihed for s xy ^ s q 


References 

[1] P. Asinari. Generalized local equilibrium in the Cascaded Lattice Boltzmann method. Phys.Rev, 
519, 2008. 

[2] A. Augier, F. Dubois, B. Graille, and P. Lallemand. On rotational invariance of lattice Boltz¬ 
mann schemes. Computers and Mathematics with Applications, 67(2):239-255, 2014. 

[3] M. K. Banda, W.-A. Yong, and A. Klar. A Stability Notion for Lattice Boltzmann Equations. 
SIAM J. Sci. Comput., 27(6):2098-2111, 2006. 

[4] P. L. Bhatnagar, E. P. Gross, and M. Krook. The Lattice Boltzmann equation: theory and 
applications. Physics Reports, (94): 145-197, 1954. 

[5] S. Chen and G. D. Doolen. Lattice Boltzmann Method for fluid flows. Annu. Rev. Fluid Mech., 
30, 329-364, 1998. 

[6] S. Dellacherie. Construction and Analysis of Lattice Boltzmann Methods Applied to ID 
Convection-Diffusion Equation. Acta. Appl. Math, 2013. 














54 


FRANQOIS DUBOIS, TONY FEVRIER, AND BENJAMIN GRAILLE 



FIGURE 37. Layout of the bounds of the L°° stability areas corre¬ 
sponding to (41c) and (41b) for the twisted scheme relative to u = V 
with an intrinsic diffusion. 


[7] P. J. Dellar. Lattice kinetic schemes for magnetohydrodynamics. J. Comput. Phys., 179(1):95- 
126, 2002. 

[8] D. d’Humieres. Generalized Lattice-Boltzmann equations. In Rarefied Gas Dynamics: Theory 
and simulation , volume 159, pages 450-458. AIAA Progress in astronomies and aeronautics, 
1992. 

[9] D. d’Humieres, Y. H. Qian, and P. Lallemand. Lattice BGK models for Navier-Stokes equation. 
EPL (Europhysics Letters), 17(6) :479, 1992. 

[10] F. Dubois. Equivalent partial differential equations of a lattice Boltzmann scheme. Computers 
and Mathematics with Applications, 55:1441-1449, 2008. 

[11] F. Dubois. Third order equivalent equation of Lattice Boltzmann scheme. Discrete Contin. 
Dyn. Syst., 23(l-2):221-248, 2009. 

[12] F. Dubois, T. Fevrier, and B. Graille. Lattice Boltzmann schemes with relative velocities. 
Accepted in Communications in Computational Physics. 

[13] F. Dubois, T. Fevrier, and B. Graille. On the stability of a relative velocity lattice Boltzmann 
scheme for compressible Navier-Stokes equations. Submitted to C.R. mecanique. 

[14] T. Fevrier. Extension et analyse des schemas de Boltzmann sur reseau : les schemas a vitesse 
relative. PhD thesis, Univ. Orsay, 2014. 

[15] M. Geier. Ab initio derivation of the cascaded Lattice Boltzmann Automaton. PhD thesis, Univ. 
Freiburg, 2006. 

[16] M. Geier, A. Greiner, and J. C. Korvink. Cascaded digital Lattice Boltzmann automata for 
high Reynolds number flow. Phys.Rev. E, 73, 066705, 2006. 

[17] I. Ginzburg, D. d’Humieres, and A. Kuzmin. Optimal Stability of Advection-Diffusion Lattice 
Boltzmann Models with Two Relaxation Times for Positive/Negative Equilibrium. J. Stat. 
Phys, 139:1090-1143, 2010. 

[18] B. Graille. Approximation of mono-dimensional hyperbolic systems : A lattice Boltzmann 
scheme as a relaxation method. J. Comp. Phys., 266(3179757):74-88, 2014. 








STABILITY OF A RELATIVE VELOCITY LATTICE BOLTZMANN SCHEME 


55 


[19] X. He, S. Chen, and R. Zhang. A lattice Boltzmann scheme for incompressible multiphase flow 
and its application of Rayleigh-Taylor instability. J. Comp. Phys., 152:642-663, 1999. 

[20] M. Junk, A. Klar, and L.-S. Luo. Asymptotic analysis of the lattice Boltzmann equation. J. 
Comput. Phys., 210:676-704, 2005. 

[21] M. Junk and W.-A. Yong. Weighted IC-stability of the lattice Boltzmann method. SIAM J. 
Num. Anal., 47(3):1651 1665, 2009. 

[22] P. Lallemand and L.-S. Luo. Theory of the Lattice Boltzmann method: Dispersion, dissipation, 
isotropy, Galilean invariance and stability. Phys.Rev., 61:6546, 2000. 

[23] P. Lascaux and R. Theodor. Analyse numerique matricielle appliquee a I’art de I’ingenieur 
(tome 1 : methodes directes). Masson, 1993. 

[24] S. Marie, D. Ricot, and P. Sagaut. Comparison between Lattice Boltzmann method and Navier- 
Stokes high order schemes for computational aeroacoustics. J. Comput. Phys., 228(4): 1056- 
1070, 2009. 

[25] M. Rheinlander. On the stability structure for lattice Boltzmann schemes. Computers and 
Mathematics with Applications, 59:2150-2167, 2010. 

[26] J. D. Sterling and S. Chen. Stability Analysis of Lattice Boltzmann methods. J. Comp. Phys., 
123:196-206, 1996. 

[27] R. F. Warming and B. J. Hyett. The modified equation approach to the stability and accuracy 
analysis of finite difference methods. J. Comp. Physics, 14:159-179, 1974. 

(Frangois Dubois) Laboratoire des Structures et des Systemes Couples, Conservatoire 

National des Arts et Metiers, Paris - departement de mathematiques, Univ Paris Sud, 

Laboratoire de mathematiques, UMR 8628, Orsay, F-91405 

E-mail address: Francois.Dubois@math.u-psud.fr 

(Tony Fevrier) LTniv Paris Sud, Laboratoire de mathematiques, UMR 8628, Orsay, F- 

91405, Orsay, F-91405 

E-mail address: Tony.Fevrier@math.u-psud.fr 

(Benjamin Graille) Univ Paris Sud, Laboratoire de mathematiques, UMR 8628, Orsay, 

F-91405, CNRS, Orsay, F-91405 

E-mail address: Benjamin.Graille@math.u-psud.fr 



